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On the Propagation of Round-Off Errors in the 
Numerical Treatment of the 


Wave Equation 
By Arnold N. Lowan 


Abstract. An upper bound of the norm of the error vector after n time steps is 
3(m + 1)(n + 2) || &* || . For the explicit scheme 5* = || &* || = 3 k 4 X 10” 
where p is the number of decimals carried in the computations. For the implicit 
scheme 6* = || 5* || is an upper bound of the errors which arise both from using 
approximations to A and A~’B in the determination of u,,, from equation (6*) 
and from rounding off the values of the products and quotients involved in the 
computation of the components of u,4; . 


Consider the numerical treatment of the differential equation of wave motion 
eu 2 eu 
1 ns SN pan 
(1) ae ° oz 


the solution of which is required to satisfy the following initial and boundary 
conditions 


(2) u(z,0) = f(z) 
(3) u(x,0) = g(x) 
(4) u(0, t) = u(a,t) = 0. 


With the differential equation (1) we will associate either of the following two dif- 
ference analogs [1] 


(5) Unear — 2a + nea = R’(urnase — 2ure + tase) 
R’ 
(6) Un ktt — 2Upn an + Una = oy (pangs — 2p ee + Une 
+ Uae — 2a + Ua et) 


where R = cAt/Az and uy, = u(hAz, kAt) with (M + 1)Az = a. 
The difference counterpart of (3) will be taken in the form 


Una — Uno 
——_— = g(hAz); 
mv) g(hAx) ; 


whence 

(7) Uni = Uno + g(hAx)At = f(hAx) + g(hAzx) At. 

The difference equations (5) and (6) may be written in the compact forms 
(5*) Un = Au — Ur 

(6*) Atiu = 4, + Burr. 


Received August 17, 1959; in revised form, December 24, 1959. 
223 











224 ARNOLD N. LOWAN 


In (5*) A is a tridiagonal matrix whose elements on the principal diagonal are 
= 2(1 — R’*) and whose elements off the principal diagonal are = R’ and u, is the 
vector whose components are the values of u(x, ¢) at time t = kAt at the lattice 
points r = hAz,h = 1,2,3--- . In (6*) A is a tridiagonal matrix whose elements 
on the principal diagonal are = 2(1 + R’) while the elements off the principal 
diagonal are = —R’ and B is a tridiagonal matrix whose elements on the principal 
oc. are = —2(1 + R’) while the elements off the principal diagonal are 


Consider first the explicit difference scheme (5*). Since both uy and u, are known, 
(5*) will yield in succession uz , u; --- . Specifically, 


Us, = Au — UW 
u; = Au, — 
(8) . 


uz, = AUn-i — Un. 


It is reasonable to assume that the components of uy are exact while those of u, , 
obtained from (7), have been rounded off to the number of decimal places to be 
carried in the computations. Let u,* denote the vector whose components are the 
rounded off values of the components of u, . It is then easily seen that we introduce 
two types of errors in the evaluation of u, . A first error is due to using u;* in lieu 
of u, . A second error is introduced as a result of rounding off of the values of the 
products involved in the expression of w+: obtained from (5) to the number of 
decimal places carried in the computations. Thus, in lieu of the exact vector uz , 
the first step in the sequence of operations (8) yields the vector u,* = Au,* — 
Uy + 5 where 5, is the error vector whose components are the round-off errors just 
discussed. Similarly, error vectors are introduced in each of the successive steps in 
the sequence of operations (8). Thus 


u.* = Au,* — uw + & 
(9) oa = Au,* — u,* + & 
un* = Aus; — ur2+5,. 
If we put ! 
(10) E, = u,* — un 
then from (8) and (9) it follows that 
(11) E, = AE, — E,n2 + 5. 


In entirely similar manner it may be shown that the counterpart of (11) for the 
implicit scheme (6*) is 

(12) E, = 4A 'E,_, + A” BE,» + 8. 

There is, however an important distinction between (11) and (12); whereas in 
(11) the components of 5, are round-off errors as above explained, in (12) the 


components of 6, are the aggregate of the errors arising both from using approxi- 
mations to A~’ and A~’B in the determination of u,,,; and the round-off errors 


he 


he 
xi- 
ors 
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introduced as a result of rounding-off the values of the products and quotients in- 
volved in the computation of the components of u, 4; . 
The error equations (11) and (12) are of the form 


(13) E, = ME, + NE,_2 + 5, . 
If in (13) we put in succession n = 2, 3, 4, --- and write 8, for E, , it may be shown 
by induction that 
(14) E, = Pai(M, N)& + Pao(M, N)& + --- 5, 
or 
(14*) E, = 2, P,(M, N)b.-» 
p=0 
where 


(15) P,(M,N) = M* + Ci_,M"°N 
+ Ci.M”“N’ + --- Ci_.M”"N’ + --- 
or 
(n/2) 


(15*) P,(M,N) = >> Cti_.M*™N’ 
s=0 


where (n/2) denotes the largest integer in n/2, where C,° = 1 and C,,” denote the 
binomial coefficient m(m — 1)(m — 2) --- (m —n+1)/nl. 

We shall prove that if M and N have the same eigenvectors, then 
(16) | P»(M, N)8,—-» || S || 8.» |] -(p + 1) 


where for any M-dimensional vector @, its norm || $ || is defined by 
1 M 

(17) Ilo ll = V@,) = 4/ 77D on” 
M h==l 


the ¢,’s being the components of 6, 
provided that the roots of the quadratic equation 


(18) a — dz — p= 0 


where the A,’s and u,’s, the eigenvalues of M and N respectively, are either 
numerically equal to or smaller than unity (if real) or have a modulus equal to or 
smaller than unity (if complex). Indeed, let 


M 
(19) t., = Daw, 


r=1 


where the w,’s are the normalized eigenvectors of the matrices M and N. From 
(19) and (14*) we get 


(p/2) M 


(20) P,(M,N)&.» = >. >, Ci. aS” ”M” *N'w, . 


s=0 r=1 
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But 
M?*N° w, = M”*y,’w, et ew, . 
whence 
(p/2) . 
P,(M, N)in-» = yaw, pT SREY 
s=0 
(21) 
= + Blp)as* Pw WwW, (say). 
It may be proved by induction that 
f) 8 —2s 8 aft ~ apt! . 
(22) Be(p) = 2) Ce-.P "u Aa Pa oe 
s=0 Vr — L2,r o=0 


where 2;,, and 22,- are the roots of the quadratic equation (18). From (22) it is 
clear that if these roots are numerically smaller than unity then 


(23) |B-(p)| <p +1; 
and furthermore 8,(p) — 0 as p — . In view of (23), (21) yields 


|| P»(M, N)8x-» || = y/ Y 1a(p)Fas PF S (+1) 4/ SD lasF 





f=) 
or 
(16) | Pp(M, N)8,-» || S (p + 1)|| Sy || - 
From (14*) the Minkowski inequality yields 
(24) || E. || s > || Po(M, N)8.-» || , 


whence, in view of (16) 





(25) IES D+ Dbl; 

and a fortiori 

5°) WE se +n = Stet?) 

where || 5* || is the largest of the sequence || 6; || , || 5: || --- || 5, |! . If 6* denotes 


an upper bound of the components of all the vectors yy , it is  endily seen that 
s* || s 3*. 


Furthermore, since 


1 M 1/2 
| E, || = -{a (Bu) 


where the £,,’s are the components of E, , it is clear that the maximum of any of 
the components is obtained by assuming that all but one of the components are = 0. 





~~ —_— he ot lL 
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Calling the maximum value of the components E,* we finally get 
(26) E,* S 3(n + 1)(n + 2)-V Me. 
The second member of (26) is an upper bound of the round-off errors for both the 
explicit analog (5) and the implicit analog (6). 
In the case of the explicit scheme (5) the matrix M of (13) is the matrix A 


appropriate to (5) while the matrix N of (13) is = —IJ where / is the M xk M 
identity matrix. The eigenvalues of A are known [2] to be 


ee) ee 
(27) ry R* cos XM +1) 
The eigenvalues of —J are clearly = —1. Thus the quadratic equation (18) be- 
comes 
2 —_ 9 =x 2 y _— = 
(28) x [2 4R* cos XM + 5 +1=0. 


It is clear that if the roots of (28) were real, one would have to be larger than 
unity, since the products of the roots is = 1. Under these conditions 8,(p) as de- 
fined in (22) would not be bounded as p — ~ and the difference scheme (5) could 
not be stable. Thus the roots of (28) must be complex, in which case the modulus 
of the roots is = 1 and 8,(p) S$ p+ 1. 

An upper bound of the round-off errors after n time steps is then given by 


E,* = 3(n + 1)(n + 2)-/Ms* 


where 6* = 3 X 4 X 10” if the computations are carried to p decimal places. 
In the case of the implicit scheme (6) matrices M and N of (13) are A and AB 
respectively where the matrices A and B appropriate to (6) have been defined 
earlier. 

It can be easily shown that the matrices A~* and A~’B have the same eigen- 


vectors, as required in the above developments [2, p. 20], and that their eigen- 
values are 


(29) » = 2/(1 + 2R’ cos’ = —] 


2 - ) : 
‘ (M 1) , ay 
Thus the quadratic equation (18) becomes 


2 


(30) — ~xz+1=0. 
1 + 2R? cos? 





— is 

2(M + 1) 

Clearly the roots of (30) must be complex. This leads to the condition 
1/{1 + 2R’ cos [rx/2(M + 1)]} <1 


which is evidently satisfied for any value of R. Thus the difference scheme (6) is 
unconditionally stable. Furthermore, and for the same reason as above, 


B(p) S pti. 


An upper bound of the round-off errors after n time steps is, therefore, once more 
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given by 
(26*) E* S 3(n + 1)(n + 2)-V MO". 


In this case, however, as previously mentioned 5* is an upper bound of the errors 
which arise both from the use of approximations to A~ and A~’B in lieu of exact 
matrices and from the process of rounding-off the values of the products and quo- 
tients involved in the evaluation of the components of u,+; . Clearly 6* depends on 
the specific scheme for solving the system of equations (6) with h = 1, 2, 3,--- M 
for the Un,b41'8. 

In order to estimate 6* for the implicit scheme we note that the counterpart of 
the typical equation (9) is 


u.* = 4A “ur, + A “Burs + & 
whence 
Au,* = 4up 4 + But. + Ad, ° 


Let R, denote the known vectors Au,* — 4us_, — Bur_,. Then A& = R, and 
therefore §, = AR, . Since the eigenvalues of A are known to be larger than 2, 
it follows that the eigenvalues of A™ are smaller than unity and therefore 


| & || = || ATRe || S| Rell. 


We conclude that 5* in equation (26*) is the largest of the norms of the n “residual 
vectors” R, = Au,* — 4ut_, — Buzz. These vectors will depend, of course, on 
the specific method of computing the u,4,’s from (6). 

A discussion of two alternative schemes for solving implicit systems of equations 
of the type (6) is contained in [3]. 
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Rotors in Polygons and Polyhedra 


By Michael Goldberg 


1. Introduction. Curves of constant width are closed curves which can be 
rotated through all orientations between two fixed parallel straight lines while re- 
maining tangent to these lines. They are also known as curves of constant breadth, 
Gleichdicke (in German), or orbiformes (in French). They have been studied by 
mathematicians beginning with Euler [1], Minkowski [3], Blaschke [8, 9], Schilling 
[11], and others to the present time. They have been applied in mechanisms to 
generate a wide variety of periodic motions. They have been used as the shapes of 
drills for drilling square and hexagonal holes. They have been produced uninten- 
tionally in certain manufacturing processes when only precise circular cylinders 
were desired. 

Two more fixed parallel lines may be added without constraining the rotation. 
In particular, the four lines may form a square. Hence, the curve may rotate while 
remaining tangent to all sides of the square. For this reason I call the curve a 
rotor in a square. 

Immediately, there arises the question whether non-circular rotors exist for 
other polygons. It has been found that they exist for all regular polygons and 
various methods of deriving them have been developed. The earliest complete de- 
velopment was published in 1909 by Meissner, a Swiss mathematician [4]. He 
derived the rotors in the n-gons and described them analytically by means of the 
polar tangential equation 


p = a + >, (a cos kd + by sin ké) 
k=1 


where p is the distance from the origin to the tangent to the curve, @ is the angle 
which the normal makes with the reference axis, the a, and } are arbitrary con- 
stants except that 


a,b. =0 for k ¥ +1 (modn). 


(For example, if n = 6, the terms which do not have to be zero are obtained for 
k = 1, 5, 7, 11, 13, 17, 19, ete.) Hence, for each polygon there is an infinity of 
different rotors. For convexity, the constants a, and b, are limited by certain in- 
equality relationships. 


2. Circular-Arc Rotors. Special cases have received special attention. In 
particular, those rotors which are bounded entirely by arcs of circles have been 
considered. Euler considered the rotor made of three equal ares, each centered on 
a vertex of an equilateral triangle. Reuleaux [2] considered rotors of an odd num- 
ber of circular ares of equal radii. They have been called Reuleaux polygons, but 
I prefer to call them Reuleaux rotors since they are not polygons. The regular 
Reuleaux rotors are rotors in a square. 





Received March 23, 1960. Retiring Presidential Address, The Philosophical Society of 
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Fig. 1.—Circular-are rotors in regular polygons. 











Fic. 2.—Regular trammel rotors in regular polygons. 


Fujiwara [14j considered rotors in a triangle. He showed that some of the 
Reuleaux rotors are also rotors in a triangle. Besides these, he derived the rotors 
of two, four and five circular ares. No such rotors of three ares exist. Of the pos- 
sible rotors of two ares, one had already been published by Reuleaux [2]. A circular- 
are rotor for the pentagon is also described by Fujiwara [14, pp. 245-246]. In 1948, 
I published a paper describing the construction of circular-are rotors for all the 
regular polygons [29]. The method is shown in Figure 1 where the centers of the 
ares are regularly distributed on two generating circles, except for n = 6 where 
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three generating circles are used. Note that for n > 4, the rotors for thé even poly- 
gons possess only one axis of symmetry (unequal generating circles), while the 
rotors in the odd polygons possess two axes of symmetry (equal generating circles). 
In these, the radii of the ares are not equal. In particular, note that the rotor in 
the pentagon is composed of two pairs of ares, the radii of one pair being twice 
the radii of the other pair. It is conjectured that, among all the rotors whose con- 
tours are made of discrete circular arcs, these have the least number of ares. 

In 1957, I published a new series of circular-are rotors which are characterized 
by higher orders of symmetry [32]. They are illustrated in Figure 2. These regular 
rotors in a regular n-gon are of two types. The upper series is made of (n — 1) 
equal segments. The lower series is made of (n + 1) equal segments. For n > 4, 
each segment is made of three or more circular ares, each are being tangent to its 
neighboring arcs except at the ends of the segments. 


3. Trammel Method of Construction. All the circular-are rotors may be ob- 
tained kinematically by a graphical method which I call a trammel construction. 
This consists of moving the rotor so that, for each portion of the motion, two points 
of the rotor trace two fixed straight lines. In Figure 3, the rotor to be generated in 
the hexagon is based on a regular pentagon shown in dotted lines. The pentagon 
is turned counter-clockwise so that the vertices B and C move along the straight 
sides of the hexagon until the next vertex D touches the hexagon. The motion is 
continued with vertices C and D moving along the sides of the hexagon. During 
these motions, the sides of the fixed hexagon will mold a rotor based on the dotted 
pentagon; that is, the rotor is the envelope of the sides of the hexagon on the plane 
of the moving pentagon. 


4. Any Rotor as the Sum of Trammel Rotors. Every rotor in a kn-gon (where 
k is any integer) is also a rotor in an n-gon. Hence, we have obtained an infinite 
number of circular-are rotors for each polygon. Furthermore, if p, = f,(@) and 
pe = f2(@) are the polar tangential equations of two trammel rotors in a given 
polygon, then their weighted mean, given analytically as 


Ps = [ufi(O) + vf2(6))/(u + v) 


where u and v are any real numbers, is the polar tangential equation of a new rotor 
in the same polygon. By extension, any rotor in a polygon can be expressed as a 





Fig. 3.—Vertices of regular (n — 1)-lobed rotor in n-gon. 
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weighted mean of a series (possibly infinite) of regular trammel rotors. This is 
similar to the Fourier series expansion used by Meissner as described in Section 1. 


5. Basic Rotors and their Construction. In Meisrner’s equation, each variable 
term in the right-hand member can serve as the entire variable part of the right- 
hand member of the equation for a rotor. These rotors, which I call basic rotors, 
resemble the regular circular-are rotors [31]. However, the curvature of the ares 
varies continuously instead of remaining constant by segments. 

Each basic rotor can be obtained kinematically by the following method. Con- 
sider a fixed circle of circumference c. Compare with Figure 4. Let another circle 
of circumference cn/(n + 1), where n is any integer, roll within the first circle 
without slipping. Then each point of the rolling circle describes a hypocycloid of 
(n + 1) cusps. Let a straight line be carried by the rolling circle. The envelope of 
the positions of this straight line will be a parallel curve of the hypocycloid. If the 
straight line is sufficiently distant from the center of the rolling circle, the envelope 
will be convex. This convex curve of (n + 1) maximal points or lobes is a basic 
rotor in a regular n-gon, as shown by the following argument. 

When the center of the rolling circle returns to its initial position, the circle 
and the carried straight line will have undergone a rotation of 2/n. As the rolling 
is repeated successively, there will be n symmetric positions of the straight line. 
Therefore, the n positions of the straight line will form a regular n-gon. The rolling 
circle can then be considered to carry the regular n-gon. As it rolls, the n envelopes 
of the n sides are the same. Therefore, each side of the n-gon keeps in contact 
with the envelope as the n-gon is rotated. Inversely, the (n + 1)-lobed rotor can 
rotate continuously within the n-gon while keeping contact with all the sides of 
the n-gon. 

A similar procedure obtains for (n — 1)-lobed rotors in a regular n-gon. 














Fia. 4.—Generation of basic rotors. 
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Fic. 5.—Conical rotors in pyramids. 


6. Rotors in Spherical Polygons. The concept of rotors has been generalized 
from the two dimensions of the plane to three and more dimensions. However, the 
intermediate case of rotors in spherical polygons has received very little attention. 
Blaschke [9] and Santalé [25] have considered ovals of constant width on the 
surface of the sphere. These are ovals which remain tangent to the two meridians 
which bound a line. 

These were generalized by the author to include rotors in all the regular spherical 
polygons. The trammel rotor method of constructing circular-are rotors in plane 
polygons is also applicable for the construction of rotors in spherical polygons. 
This time, however, the generated arcs are not plane and, therefore, they are not 
circular [30]. 

Also, the spherical rotors corresponding to the basic plane rotors can be gen- 
erated by the rolling circle method. In this case, the rolling circle carries an arc 
of a great circle instead of a straight line. (See Figure 4.) The symmetric positions 
of this great circle make a regular spherical polygon [31]. 

Every rotor in a plane polygon has its counterpart as a rotor in a spherical 
polygon. However, there are more types of spherical rotors for two reasons. The 
most obvious reason is that regular spherical polygons of the same number of sides 
are not similar; their shape depends upon their size. Therefore, their rotors are 
correspondingly different. But the more surprising difference between plane and 
spherical rotors is the fact that spherical ovals of constant width, which are tangent 
to two ares of great circles, are distinct from rotors in a spherical quadrilateral. In 
the plane, these two sets of rotors are identical. 

Models of spherical surfaces are difficult to construct and to demonstrate. 
However, every spherical polygon can be converted into a pyramid by passing 
planes through the great circles by which it is bounded. The rotor in the spherical 
polygon is replaced by a non-circular cone which is now a rotor in a pyramid. 
Models of several conical trammel rotors are shown in Figure 5. 


7. Rotors Tangent to n Circles. The trammel method and the rolling circle 
method may be applied in the derivation of other types of rotors. For example, 
n equal circles may repla:e the n straight lines of a plane regular n-gon [33]. Among 
the interesting possibilities is a series of rotors which approximate regular polygons. 
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Fig. 6.—Rotors tangent to n fixed circles. 


Several are shown in Figure 6. Also, the methods are just as applicable on the 
surface of a sphere. 

A further generalization is the use of a set of n symmetric curves instead of n 
circles. No special interesting cases are available at present. 


8. Rotors in Regular Polyhedra. The first extensive study of surfaces of con- 
stant width was made by Minkowski [3]. He devised general geometric methods 
of deriving many of them. Many beautiful theorems concerning them were ob- 
tained. The most obvious method of generating one is by revolving a symmetric 
oval of constant width about its axis of symmetry. A more unusual surface is based 
on a regular tetrahedron [6]. Each vertex serves as a center of a spherical surface 
passing through the other vertices. However, one edge of each pair of opposite 
edges must be chamfered to a portion of a toroidal surface. The foregoing shapes 
are shown as A and B of Figure 7. 

The surface of constant width is a rotor in a cube. The general investigation of 
rotors for all the regular polyhedra was first considered by Meissner. In a very 
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elegant paper [7], he showed that non-spherical rotors exist for the regular tetra- 
hedron, the cube and the regular octahedron, and that they do not exist for the 
regular dodecahedron (12 faces) and the regular icosahedron (20 faces). Just as 
Meissner used a Fourier series development for rotors in plane polygons, he used 
a spherical harmonic development for studying rotors in regular polyhedra. He 


showed that the rotors in the cube may be described by the polar tangential equa- 
tion 


p(6,¢) =~aot+¥it+¥3;+¥3+¥74+-::-; 


the rotors in the tetrahedron may be written as 


p(6,¢) = ao + Yi,+ ¥2+ Ys; 


and the rotors in the octahedron may be written as 
p(6,¢) = a+ ¥i4+ Ys, 


where Y; (a function of @ and ¢) is the spherical surface harmonic of the ith degree. 

The number of arbitrary parameters in the equation for the rotor in the cube 
is infinite. The number of parameters for the rotor in the tetrahedron is eleven, 
while the number of parameters for the rotor in the octahedron is eight. Several 
rotors for the tetrahedron are shown in Figure 8. A rotor for the octahedron is 
shown on the extreme right of Figure 7. 

Before leaving three dimensions, it may be fitting to consider other possibilities. 
It is conceivable that non-regular polyhedra may have non-spherical rotors. In 
particular, it is not known whether open-ended regular prisms can have non- 





Fie. 7.—Rotors in cube. 
A. Rotor based on tetrahedron; B. Rotor of revolution; C. Three-lobed rotor; D. Five- 
lobed rotor; E. Rotor in octahedron. 





Fig. 8.—Rotors in regular tetrahedron. 
A. Prolate rotor of revolution; B. Oblate rotor of revolution; C. Triaxial rotor 
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spherical rotors besides the weil-known rotors in a rhombic prism which are surfaces 
of constant width. 


9. Rotors in Higher Dimensions. Surfaces of constant width have been gen- 
eralized to rotors in higher dimensions. Santalé [26] derived relations between 
measures of their boundaries and their contents as generalizations of relations 
derived by Minkowski [3, pp. 215-220]. Rotors in the simplex (the analogue of 
the triangle and the tetrahedron) also exist. Simple examples are the bodies of 
revolution given by the following polar tangential equations: 


p = a+ bcos2¢, Rotor in simplex, 
p = a+ bcos3¢, Rotor in hypercube. 


Non-spherical rotors for the analogue of the octahedron have not been mentioned 
in the literature. 


10. Polygon Rotors in Ovals. If a rotor in a regular polygon is held fixed while 
the polygon is rotated about it, all the vertices of the polygon trace the same 
curve. Therefore, if this curve is fixed, the regular polygon can be rotated within 
it while all the vertices lie on the curve. Hence, for each rotor in a polygon, we 
have a curve within which the polygon can be rotated. See Figure 9. An application 
of a triangular rotor in an oval is a recent design, by the German engineer Felix 
Wankel, of a non-reciprocating internal combustion engine. Test models of this 
engine have been made for the Curtiss-Wright Corporation. 

This relation also applies for rotors in spherical polygons from which we obtain 
pyramidal rotors in non-circular cones. 

However, in three dimensions, this relation does not apply. If a non-spherical 
rotor in a polyhedron is held fixed while the polyhedron is rotated, the vertices 
of the polyhedron do not lie on a surface. Instead, the positions of the vertices 


= 

















eo 


Fie. 9.—(a) Triangle rotor in an oval; (b) Square rotor in an oval. 
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Fic. 10.—Double-contact cam. 


fill a volume bounded by two closed surfaces. The ellipsoid is the only known 
surface possessing the property that the vertices of all the circumscribing rectangu- 
lar parallelepipeds lie on a sphere. It has been conjectured that every other body 
generates a volume in this manner and not a surface [10]. 


11. Mechanisms Related to Rotors. The use of drills in the shape of rotors in 
polygons has been extended from the drilling of even polygons (the square and 
the hexagon) to the drilling of odd polygons (the triangle and the pentagon). 
Several other problems, closely related to rotors in polygons, have been investi- 
gated. One is the determination of the non-circular shapes of pivoted rotors which 
remain in contact with two straight arms of a pivoted rocker arm. See Figure 10. 
Again, the only admissible angles between the arms are rational fractions of a 
circle as are the angles of a regular polygon. Examples of the ovals are described by 


pi(@) = cos p2(8) 


where p2(0) = mx/2n + k sin n@, which is a rotor in a polygon. A geometer would 
describe these ovals as curves whose isoptic curves are circles [28]. They are the 
basis of a patent issued for a series of double-contact cam mechanisms [36]. 

Another related mechanism is the intermittent rotor [34]. This makes contact 
with a series of fixed elements but not always with all the elements. In the example 
shown in Figure 11, the rotor is restrained in its motion by contact with three of 
the four fixed points until all four of the fixed points are touched. The motion may 
then be continued with another set of three fixed points as constraints. 

The presentation of new and unfamiliar basic mechanisms in this paper shows 
that the ancient science of mechanisms is far from exhausted. When the engineer 
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Fig. 11.—Intermittent rotor. 


and designer become more familiar with them, it is expected that more applications 
will be made. 
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Numerical Quadrature Over a Rectangular 
Domain in Two or More Dimensions 


Part 3. Quadrature of a Harmonic Integrand 


By J. C. P. Miller 


1. Introduction. In Part 1 [1], §5, formula (B), §7, formula (B’), and in §9; 
also in Part 2 [2] in several places, we have seen how the error term is very much 
reduced if the integrand f(z, y) is a harmonic function, that is, if V’f = 0. In this 
note we pursue further this special case, in which especially high accuracy is at- 
tainable with few points. 

It may not be often that the integrand will have this special form, but it seems 
worthwhile to develop a few of the interesting formulas. We start by obtaining 
expansions for n variables, and more extensive ones for two variables, and then 
obtain and consider special quadrature formulas. 


2. Expansions. As in Part 2 [2] §2, we develop f(x , t2,--- , t,) as a Taylor 
series in even powers of each of the variables x, . Then, using Vf. = 0 whenever it 
is applicable, we obtain 


la 


























= [/(2h)" = ny ( LY sa, Xn, °*+, tn) da, da. --- dx, 
= fet soya h Boy 4 5ier+ Pe) h 
(2.1) < + a (13 ov 7 ’) f. 
+Gastos ton + ae 
| hs 


where, as before, extended, 


a) he-Doe, we Doe, hh - 5S 


02,70x,? 02,702,202? 02,202,202 202.2 


etc., the summations extending over all possible combinations of r, s, t, - - 
no two equal. 


Labelling the symmetrical sets of points as in Part 2, we have likewise the ex- 
pansions for sums of values of f over the sets 


- with 





Received December 30, 1959. This work was supported in part by the Office of Naval Re- 
search. 


240 


a = we 


ith 


ex- 


Re- 


NUMERICAL QUADRATURE OVER A RECTANGULAR DOMAIN 241 


(231) 0 fo 








4hé 4 n°, 6 8 8 
(232) ala) 2nfo— AS wif, + TO otf, + AO (at — 26" Wf 
1 0 
ai Sd a (p's ® °”)f, 
12 12 
=a ae (2p" — 35” 


> 6D‘e° + 68") fo + pee) 


(233) A(b) n(n — Ifo — 2 (n — 4) 0%, + “E? (n — 16)shfe 


8,8 
+ add {(n + 6)D* — 2(n — 64)Q"}fo 
10; 10 
‘ cD {(n — 4) p's" — (n — 256)0"}fo 
12; 12 
- er {2(n — 34)” — 3(n + 362)3” 


— 6(n — 728)D'e° + 6(n — 1024)8"}f. + --- 
(2.34) y(e,d) 4n(n — 1)fp — = {(n — 1)(c' + d*) — 6c'd"} Df, 
6 
+ sl {(n — 1)(c° + d*) — 15e'd*(c* + d*)}5°fo 
sh® 8 ‘ 23727 4 =0,.4 s 
+3 [{(m — 1)(e + d*) — 28c'd*(c* + d*) + 70c'd'}D 


— 2{(n — 1)(e + a) — 28cd*(c' + d*) — 70c'd'}e'lfo 








eee 
4 DRE (nt — 33n + 122)stf 
+ oi {(n — 2)(n + 13)D° — 2(n* — 129n + 1094)Q*}fo 
ries 


We recall that 0 is the origin, or centre of the square, a(a) includes all points with 
one coordinate --ah and the rest zero, 8(b) has two coordinates each independently 
+bh and the rest zero, y(c, d) has one coordinate +-ch, another +dh and the rest 


zero, and finally «(e) has three coordinates each independently --eh with the rest 
zero. 
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3. Expansions over a Square. Such expansions are simpler since 3°fy, Q°fo etc., 
are absent. They can be obtained by analysis with the detached operators—in 


particular 2; we proceed to obtain expansions with general terms. 


If F(z) = u + w is a function of a complex variable z = x + iy then both 
u and v are harmonic functions satisfying D,@ + D,’¢ = 0. Likewise, if u is a 
harmonic function, it can be shown that v exists such that u + iv is a function of 


a complex variable. We then have 
D,F = iF’ = iD.F 











and 
(3.1) D.D, = D = iD? = —iD,. 
In order to develop expansions we therefore substitute 
(3.2) D, = i>" D, = ig, 
Consider, firstly 
h h 
(33) J = (2h)~ [, [f y) dx dy = a (eh?= — es) (ety — eo *r) fy. 
z u 
The operator is 
sinh hD, sinh hD, _ sinh i ‘“h® sinh i'"hD 
h?D. D, " h?D? 
(3.4) et 1 cosh (a? + i-*)hD — cosh (2°? — i" )AD 
sig 2 h?D? 
1 cosh /2h D — cos /2hD 
2 hp? 


whence 


2: ai ‘ad 


(3.5) i= (2 +i 4 2h os, + Po +) he, 


10! (4r + 2)! 
Likewise 


» f(a, Ya) 


(es + g Pe 7 ane + e ”) f, 
2 (cosh ahD, + cosh ahD,)fy 
2 (cosh i "ahD + cosh i*ahD)f, 


(3.6) ah ah 
= 4 cosh < 75 D 008 75 Dfe 
ah .., avh* ah” 
=s[1-ar4 Se +(-) Gym" 
and 
{ > flag, ys) = 4 cosh bhD, cosh bhD, fo 
B(b) 
= 4 cosh ¢ "bh®D cosh 7"bhD fy 
(3.7) 4 = 2 (cosh bhv/2D + cos bhv/2D)fy 





=af1 2’ “ x 2 vn 2 oh” 





(4r)! 


me, 
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eee b= = a” +--+ |. 
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We shall not use all the expansions given above in the present rote, but it 
seems useful to set out the collected results for future use. 


4. Lattice-point Formulas over a Square. We consider first formulas in two 
variables, and start with 9 points, putting a = 6 = 1 and using the sets 0, a(1), 
B(1). We write 
(4.1) J = 1/4W = Aof(0,0) + Do Acf(ta, Ya) + dD Aaf( zs , yp) 
using (1a, Ya) etc. as typical sets of coordinates. 

Using (3.5) to (3.7), we equate coefficients of Df), r = 0(1)2. This gives 


Ap +4Aat+ 4Ag=1 
(4.2) —4A,+ 16 Ay = 75 
+4A,+ 664A; = }$ 


» 12 
with correction term C= -(-44. + 256A, — *) a Dfo. 
We obtain the formula 





| 


| 7 —32 v | + 900 
; ; 1952 h” | 1 
(43) | —32 1000 —32 | with main correction term — 1365 121 Py. 
7 —32 7 | 








This formula is remarkably good. With the example of Part I, we have, writing 
J’ = WJ 


1 1.2 1.2 1 
J’ = il [ sin z sinh y dz dy ri (1 — cos 1.2) (cosh 1.2 — 1) 
0 


0.12922 70590 73675 11602 


Formula (4.3) gives 


J’ = 0.12922 70590 72834 11029 
with EF = —0.0°841 00573 and C = +0.0°841 01633. 


5. Five-point Formulas. The high precision of (4.3) suggests that formulas of 
lesser precision, with fewer points, may be useful. We use the first two of (4.2) 
and take one of Ao, Aa, Ag to be zero. 


(i) Ao = O gives an eight-point formula with relatively poor precision. 


19 56 19| + 300 


’ : . h* 
(5.1) 56 0 56) with main correction term — 40 Fy D'fo . 


19 56 19 | 





11 0 1/+60 
| | 29 8 
(5.2) 0 56 0 | with main correction term — : = D'fo . 
a aa 
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(iii) As = 0 gives 





on 0|+ 15 


8 
—-1 19 —1 | with main correction term + ap D*fo. 


Sodas 0 | 


(5.3) 


We observe that (5.2) and (5.3) combine in the proportions 5:7°; to give 
(4.3), though without an error estimate! Likewise $ X (5.2) — $ X (5.3) gives (B) 
112 h® 


of Note I, and an estimate for the correction, namely — —— rT 
o 3 


D*fo when f(z, y) 


is harmonic. 
Another combination, that of (5.2) and (5.3) in equal proportions, gives a 
small correction term: 


a —-—; 


1 -4 1/+ 120 
| ; , 4h’. 
(5.4) —4 132 —4| with main correction term — = 51 Dfo. 
og: 
1 —4 1 | 


Again 4 X (5.2) — 3 X (5.3) gives small multipliers: 


1 3 1/+15 

3 -1 3| 

1 3 1 | 
Evidently (4.3) is most precise, but simultaneous use of (5.2) and (5.3) gives 

an idea of the precision attained, and readily yields the better result if desired. 


Formula (5.5) might be helpful with desk computing, but (5.1) has little to recom- 
mend it. 


Numerical results for some of the formulas using the example of §4 are as follows: 


9 
with main correction term — — 


o~ 
o 
tn 

~~" 














Formula Result J’ 10° X E 10° x C 











= — —_—n — ; " a = 
(5.1) | 0.12922 72986 +2395 —2396 
(5.2) 0.12922 70974 +383 —383 
(5.3) 0.12922 70255 —336 +335 
1(B) | 0.12922 71932 +1341 —1342 


(5.4) 0.12922 70615 +2 | —24 





6. General n; 2n? + 1 Points. We consider now the n-dimensional case, n = 3, 
using lattice points 0, a(1), 8(1). In this case the term in 3°fy is relevant, and the 
©*fo term will appear in the error, except when n = 3. 

We equate coefficients of fo , D'fo , 3°fo in the expansions resulting from use of 


(2.1), (2.31)—(2.33) in (4.1). We obtain 
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Ao + 2n Ag + 2n(m — 1) Ag =1 
(6.1) — 4A,— 8(n—4) As = #5 
+ 6A, + 12(n — 16) Ag = $$ 

while 


8 
C= -{ 44. + 8(n + 6)4, — 1K oy 


92 | a’ 
cs { 8A, + 16(n — 64)Ag + +} gi Se - 




















These yield 
—61n’ + 931n + 3780 61n — 496 61 
ao Ana. dom — 
(62). Me 3780 . 3780 ts 7560 
with 
, 1198h° 5 | 3619h" 
C= om gi Dt 315 sie!" 
In particular 
12048 626 61 
3. = =o kena ws ~aom 
GH st AGE 64 7500 * — — 7560 
13056 504 61 
J = he a ane Ac = 7- >> Z G8 @ ese 
— Hes 75 4 70 7560 
, i. _ 13820 a 
(62) s-5 Aa (4 ae ee 
_ 14340 eS 
(6.36) n=6 Ay = 7560 A, = 7560 Ags = 7560 
As a numerical illustration for n = 3 we consider 
1 Va del 5 
=-l=-; 0S — osh — z dx dy dz 
J 3! 51 Lf, 0852 08 y cosh 3 2 dz dy 
= ba sin 4 sin 1 sinh 2 = 0.9800827. 
15 4 4 


Formula (6.33) gives J = 0.9799734 with E = —0.0000109 and C = +0.0000110. 
This result is less spectacular than that of §4, for these reasons: 
i) In $4, h = 0.6, here h = 1, and the correction term in (4.3) contains a high 
power of h. 

ii) The correction term in (6.2) is of order h*, that in (4.3) is of order h” 

iii) The higher the number of dimensions, the more individual terms there are in 
D*f, Df, etc. In (4.3) there is only one term in D"f, in (6.33) there are 9 in 
D'f. 

iv) The effect of larger interval Ah is enhanced by the use of the factor }, which 
exceeds unity, in cosh } z; this is only partially balanced by the factor cos } z. 

In spite of these points, the formula (6.2) seems a good one. 


7. Quadrature over a Square; Specially Chosen Points. Since the expansions 
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of §3 involve only cross-differences D‘fy , it appears likely that use of sets of diag- 
onal points 8 will be more profitable than attempts to use sets a. It turns out that 
sets 0, a(a), 8(b) and 0, a(a) both fail to give real values of a if maximum preci- 
sion is sought. On the other hand, we can get several formulas making use of any 
number of sets 6(b,), p = 0(1)r, both with and without the point 0. 

We start first with r sets 8(b,), without the point 0. We have to find the 2r 
constants Ag, , b, satisfying the equations 

1 


. 4(s—1) - an - 
(7.1) > 4Ag, by @e— Gens =O 8 = 1M) 


obtained by substitution of (3.5) to (3.7) in 
(7.2) J = > As,f(+b,h, +b,h) 


and equating the coefficients of the first 2r coefficients of D*. Sundry powers of 4 
have been cancelled. 
By familiar arguments, the b,' are roots of the equation 





kis x Te 
Co Ci C2 iS a 

(7.3) C; C2 C3 tl Cosa = 0. 
Ge Cra C,-2 pe Coa 


These are the orthogonal polynomials for the weight function w(x) = 3(a** — 2”) 
and range 0 < x < 1. The first two are 


(152 -—-1=0 
(7.4) 
\s1922 — 4384 + 11 = 0. 
The main correction term is obtained from the next power of * and yields 
r ‘ on D*” 
5 C=(|C, — 4Ag bs } ——— ho. 
si (Co — 3 44,00) 


If the point 0 is included, our equations (7.1) are replaced by 


| do + 4>° As, = 1 
p=l 





(7.6) i 
| 2d, 4Ag, by = mes nae 5 C,, $s =1(1)2r 
and the b,' are roots of the equation 
1 x x oo 
1 Ci C. C3 ee OPES 
C. Cs Cy +++ Crel/=0 


Y 
Cri Cr+2 C43 te Cy 
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which are the orthogonal polynomials for the weight function w(x) = 4(2"* — 2”) 
for the range 0 < x S 1. The first two are 


3x —-1=0 
(7.8) 
17017x° — 13650z + 1745 = 0. 


The main correction term is this time 
ofrte ote 


wt aay : 8r+4 ime. 
(7.9) C= (Cons 2d, 4Ag, by ar Di 


In each case the coefficients As, may be computed by standard methods. 


a” fo : 


8. Formulas for r = 1. These have 4 and 5 points respectively 
64 h*p* 


= — = 54 = — —_— 
(8.1) asin Tiree: C= 35 81 
(8.2) Ay = 4 4e = 1 b=a Zz C= 2816 nD" 
- a ne 12285 12! “°° 


Written out in full: 
(8.3) J=H = HS, 15) + f(-15" 4, 15) 
+ fs, -15*") + f(-15""™, -15"™)} 
4f(0, 0) + Hof f(3, 3) + f(—-3™, 3) 
+ f(s", —3-) + f(-—3™, —3™)}. 


(8.4) J 


Il 
a 
— 

Il 


As a numerical test use 


1 1 
J= il = j [ | cos x cosh y dx dy = sin 1 sinh 1 = 0.98889 77057 62865. 
1 #1 


Formula (8.3) gives 0.98889 06525 with 
E = —0.00000 70533 and C = +0.00000 70547 
and formula (8.4) gives 0.98889 77062 41358 with 
E = +0.0°4 78493 and C = —0.0°4 78543. 


9. Formulas for r = 2. These have 8 and 9 points respectively 


b; = 0.40316 26030 59346 89754 Ag, = 0.22912 30654 28169 97222 

(9.1) 
be = 0.84439 75319 23478 74713 Ag, = 0.02087 69345 71830 02778 
; ; . 54592 256 
with main correction term 57014685 4 16! 
bo = Ap = 0.69521 80834 12925 81989 
(9.2) b, = 0.63205 02078 18796 99524 Ags, = 0.06686 42185 46105 38162 
be = 0.89531 63791 24106 97730 Ag, = 0.00993 12606 00663 16340 


h*a"*f, 
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16832 x 1024 haf, 


with main correction term 78975897 201 


For 
1 1 1 
J=;/ [ cos x cosh y dz dy 
4Ji ti 


formula (9.1) gives 0.98889 77057 62853 38396 with 
E = —0.0° 1171243 and C = +0.0"1171555 
while formula (9.2) gives 0.98889 77057 62865 09647 with 
E = +0.0°9 and C = —0.0"90. 
With formula (9.2) we find approximately 0.82447 37090 77903 16756 for 


2 2 
a [ [ cos x cosh y dz dy = sin 2 sinh 2 = 0.82447 37090 77809 15433 
2 2 


with EF = +0.0°9401323 and C = —0.0%9406250. 
These formulae clearly have high precision, even with considerable values of h. 


10. Quadrature over a Cube; Specially Chosen Points. The search for such 
formulas is more difficult in 3 or more dimensions. It seems that one or more extra 
available constants are needed in order to obtain real points. We shall not pursue 
this, but give one simple formula for three dimensions. 

We find nothing convenient by use of points a(a), with or without 0; likewise 
0 with 6(b) fails to give real points. We can, however, use 12 points 6(b) alone. 
We have then to satisfy 


(Qn(n — 1)Ag = 1 
(10.1) 
\8b°(4 — n)Ag = v's, where n = 3. 


This yields b = (2/5)"* = 0.79527 07287 67051 Ag = 1/12 


with main correction term 


x 6 16\ h° 


a 3°fy = 0.005626h'S'f,. 

With the example of §6, with integrand cos $x cos y cosh $z (10.1) gives 
J = 0.97519 with E = —0.00489 and C = +0.00494. 

The only formula found that allows for the term 3°f) and has an error of order 
h* is (6.33), which needs 19 points. It is evident that further search is needed. 
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A Method for Calculating Solutions of Parabolic 
Equations with a Free Boundary 


By Milton E. Rose 


1. Introduction. The description of phenomena involving phase transitions 
leads to a class of parabolic partial differential equations with a free interior bound- 
ary which marks the interface separating the phases. The simplest nonlinear prob- 
lem of this type was first treated by J. Stefan in 1891 and considerable attention 
has been focused on such problems in recent years [4], [6], [9]. 

A typical mathematical formulation of this type of problem is the following: 

Prosiem I: Determine functions x(t) and u(z, t) satisfying 


(a) tUse(z, t) = u,(z, t), 0 < z < x(t), t>o 


(b) u.(0,¢) = —1, t>0 

(1) (c) #(t) = —u,(x(t), t), t>0o 
(d) u(z, t) = 0, = #6), t20 
(e) 2(0) = 0. 


With suitable choice of units this formulation describes, for example, the recrystal- 
lization due to uniform heating at one end of a semi-infinite substance which is 
initially at its critical temperature [3]. Equations (1c), (1d), and (le) serve to de- 
scribe the interface separating the phases of crystalline state. Effective numerical 
procedures and proofs of their validity have been described by Douglas and Gallie 
{2}, and Trench [8]. 

The object of this paper is to suggest a computational approach to general 
problems of this type which has the feature that the path of the interface is not 
regarded as an explicitly imposed interior boundary condition. The method is 
closely related to a proposal by P. Lax [7] for calculating weak solutions of hyper- 
bolic equations in conservation form as the limits of suitable finite difference equa- 
tions. Several problems involving shocks and contact discontinuities have been 
successfully treated by Lax from this point of view. The present investigation was 
motivated by the question of whether or not phase-change interfaces were sus- 
ceptible to a similar treatment, particularly in view of the fact that such interfaces 
also arise as discontinuity surfaces in the general treatment of the equations of 
hydrodynamics (Keller [5)). 

In the first part of this paper a formulation of certain phase transition phenomena 
is described (Problem II) and the concept of their weak solutions is introduced. 
Heuristic arguments are given serving to identify a weak solution of Problem II 
with the solution of Problem I. The results of several calculations are cited to lend 
support to a conjecture that the solution of certain natural finite difference equa- 
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tions related to Problem II converge to a relevant (essentially unique) weak solu- 
tion of the problem. In addition, a random walk process suggested by the difference 
equation is mentioned which tends to illuminate certain aspects of phase transition 
processes and suggests an additional computational approach to such problems. 

The author wishes to express his appreciation to P. Lax and E. Isaacson for 
their interest and discussions concerning this paper and to N. Metropolis for dis- 
cussions of the random walk formulation. 


2. Weak Solutions of Phase Transition Processes. We begin by deriving equa- 
tions governing a particular example of a process involving a phase transition.* 
Consider [3] a bar of uniform cross-section of a substance which undergoes a change 
in crystalline form with negligibl. change in density involving a latent heat of 
crystallization. Let e denote the specific internal energy, q the heat flux, 7 the tem- 
perature and p the density. The equation 


(2) (pe): + ge = 0 


expresses the conservation of energy of the process (we use the notation v, = 
dv/dzx, etc.). 

We assume Fourier’s law, 
(3) q= —kT, 


where k is the coefficient of thermal conductivity. Finally recalling our assumption 
that the density is constant, we assume an equation of state given by 


(4) T = T(e). 
In a region in which this transformation is one to one we obtain 
T. = aT es 
with a = k/per. This is the usual equation of heat conduction. However for 


processes in which a latent heat of recrystallization occurs the transformation (4) 


is no longer one to one; the following relationship is typical of such processes (see 
Figure 1). 


(4’) T = y(e) + To 
with 
y ‘(e — @), e< &% 
y(e) = 40, @oSeSsat+H 
y*-[e — (H + &)], e> e+ H. 


Here H is the latent heat of recrystallization and y and y~ are non-negative 
quantities related to the specific heat of the material. Without loss of generality we 
may assume é = 0 and J) = 0 

If we adjoin to equations (2), (3) and (4’) appropriate initial and boundary 
conditions it seems reasonable to suppose our formulation describes a unique and 
well-posed physical problem. Summarizing, we formulate 


* Our discussion is easily generalized to several space dimensions. 
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Fic. 1.—Equation of state for a recrystallization process. 


ProsLeEM II: Find solutions e(z, t), g(x, t) and T(z, t) of the equations 


(2) (pe): + ge = 0, z>0,t>0, 

(3) q+kT, = 0, z>0,t>0 
(ve, e<0 

(4”) T = 40, O<sesH 
lye a H), e> H 

satisfying the initial and boundary conditions 

(5) e(x,0) = e(z), z>0 

(6) q(0, t) = g(t), t>0o0 


Because of the mild nonlinearity introduced by equation (4”) our formulation 
of Problem II may not allow a genuine solution, i.e., a solution which is continuous 
with continuous derivatives in the quarter plane D [x > 0, t > 0] for a compact 
set of initial and boundary data. The physical situation suggests that the class of 
solutions be enlarged to allow jump discontinuities along certain smooth curves 
C in D. 

We do this by introducing in the usual manner the concept of a weak solution 
of Problem II. Consider a space @ of pairs of smooth testing functions (¢, ¥) with 
compact support, i.e., each (y, y) in ® vanishes identically outside some compact 
bounded region which lies in the half space ¢ > 0 and does not intersect the bound- 
ary x = 0. We multiply equation (2) by ¢ and equation (3) by ¥, integrate and 
then integrate by parts to obtain 


(7) [[ Get+eca)dzat+ f o(2,0)eo(2) dz = 0 


and 


(8) [| (qy — kTy.) dxdt = 0 
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recalling that e(z,0) = eo(xz). We restrict (e, g, JT’) to the class of functions satisfy- 
ing the boundary condition (6) and the equation of state (4”) and call (e, g, T) a 
weak solution of Problem II if (7) and (8) hold for every set of test vectors (y, w) 
in ®, 

A genuine solution is, clearly, a weak solution. On the other hand, two genuine 
solutions whose domains of definition are separated by a smooth curve C given by 
zx = x(t) will constitute a weak solution if and only if the following equations 
(obtained by applying Green’s theorem) hold: 


[pel = [a] 
and 
[T|z = 0 


where [v] denotes the jump in the quantity v across C and ¢ =, dx/dt. Using (3), 
(4”) and the assumption that p is constant we obtain 


(9) pHz = —(kT) 
and 
(10) [T] = 0 


which express the conditions of balance of heat flux and the continuity of tempera- 
ture across the interface. These conditions correspond to the conditions expressed 
in equations (1c) and (1d) of Problem I; they result in Problem II as a requirement 
for a weak solution. 

I conjecture that a weak solution of Problem II exists which is composed of 
piecewise genuine solutions; furthermore it is the only such solution which is con- 
tinuously dependent on the initial conditions with respect to some appropriate 
norm. If this conjecture can be proved correct, the identification of this solution 
with the solution of the Stefan problem, Problem I, would follow directly through 
proper choice of the initial and boundary conditions. Probably the solution of Prob- 
lem II can be obtained as the limit of a certain set of finite difference equations 
which are described below. More precisely, 


ConsEecTurE A: There is a unique weak solution (e, g, 7’) of Problem II com- 
posed of genuine solutions in each of the regions separated by the interface curve 
x = x(t); furthermore, the solution uw of Problem I and the solution (e, g, 7) of 
Problem II (with suitable initial and boundary conditions) as identified through 
the relationship 7'(e) = wu are identical. 

ConsEcTuRE B: (e, g, 7’) may be obtained as the strong limit of the sequence 
of solutions of the finite difference equations described below (Problem II,). 


I wish to give certain plausibility arguments which, together with the results 
of certain calculations described later, tend to support these conjectures. In doing 
so I lean heavily on the already cited paper by Lax [7]. 

The following describes the simplest finite difference approximant to Problem 
II: at a point (7, n), 7 = 1, 2,--- ,n = 1, 2,---, of a mesh A of size (Az, At) 
replace e; by (e;"** — e;")/At, ge by (qi” — gi-s)/Az, and T, by (T?, — T,")/Az. 
This leads to 
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Prosiem II, : Find solutions e;”, ¢g;", and 7;" of the equatioris 


(24) en = e;" = u( qi” cag qi) 
(34) qi” = —k( Thy > a. ? 
(y~e;", e,” <0 
(44 T;* == 0, 0 <= e;”" s H 
y*(e:" — H), ce" >H 
satisfying the initial and boundary conditions 
(5a) e, = €o(tAx) = Zoi 
(64) go” = g(ndt) = g” 


Here, u = At/pAx and k = k/Az. 
Let us suppose initial and boundary conditions for which ¢; = 0 and g” = 0. 


It is not difficult to show then that e;" = 0 and the above equations can be simpli- 
fied. Letting 


»_f0 if O<e*<H 
Pi “V5 if e; >H 


and eliminating q;" and T;" in (24), (34) we obtain 


(ef — H) = piu(etn — H) + piu(etu — A) 


(11) 

— (1 -> 2 pi") (e:" —_ H) 
forn = 0,7 = 2, and 
(12) (e?** — H) = po"(eo” — H) + (1 — pr") (eX” — A) + why” 


forn = 0. 


Very likely the stability and convergence of this system is insured by requiring 
that the coefficients (1 — 2 p;") be non-negative. This leads to the condition 





At p 
(13) at © Bet 


lA 


From the definition of p;” it is easily seen from (11) that, if e74:, e;", and e?_, 
are each in the interval [0, H], then e?7** = e,;". Suppose, for example, that ini- 
tially e; = 0. We may conclude that the width of the transition region between 
phases, i.e., the region in which e;" varies from its initial value 0 to a value greater than 
H, is 2Az for any line t, = constant. This property provides a striking contrast to 
the spreading of shock zones in Lax’s treatment of hydrodynamic problems and 
suggests that phase transitions can be determined accurately by this method. 

It is of interest to consider the interpretation of equation (11) as a random 
walk process. Consider wells each capable of holding H balls at points Py , Pi, --- 
situated on a line. A ball in a well which is not yet filled is required to remain in the 
well (p;” = 0); a ball at a point at which the well is already filled moves to the 
right or left with equal probability p;" > 0 or remains at this point with probability 
(1 — 2 p,”). 

The boundary condition qo" describes the rate at which balls enter the process 
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at the left and the capacity of each well measures the latent heat of the process. 
In addition to their theoretical interest, random walk processes of this type suggest 
feasible computational methods for more general, but related, problems. 


3. Numerical Computations. In this section the results of several calculations 
based on the difference scheme outlined in Problem II are described for the following 
particular problem: 


é+q.=0 
q+T.=0 
rai? Os e<1 
de, e<0 

In each case the initial and boundary data were 

e(z,0) = 1, x>0 
and 

q(0,¢) = —1, t> 0. 


The interface curve x(t) was determined as the distance from the origin to the 
first meshpoint on the line ¢, = constant for which 7’ = 0. For comparison (c.f. 
eq. (9)) 
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Fia. 2.—Interface for Ax = .02, \ = .5. 
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Fig. 3.—Dependence of relative error on X. 


TABLE 1 
Comparison of values t(x, Ax), (At/Az? = 4) 





Ra 


| u(x, .1) (x, .02) 








j 0.110 a4 0.108 


0:1 | 

0.2 | 0.225 0.232 
0.3 | 0.365 | 0.370 
0.4 | 0.515 | 0.519 
0.5 0.675 0.679 
0.6 | 0.845 in 
0.7 | 1.025 | sail 
0.8 


| 1.220 —* 


* Values not calculated. 








was integrated numerically using a backward difference for T, at x(t) and the per- 
cent relative error 


a = 100 [2= 200] 
x(t) 
calculated. 

Figure 2 shows the interface for the case where Ar = .02 and At/Az’ = .5. 
Comparison to the case in which Az = .1 and At/Az* = .5 may be made from 
Table 1 which compares i(z; Ax) for these two cases. 

In Figure 3 the dependence of the relative error d on the value \ = At/Az’ is 
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given. In each case Ar = .02. As expected, instability was observed for values of 
\ > 1, although this fact has not been indicated in the figure. 


4. Conclusions. The numerical evidence presented here, although somewhat 
limited, seems to support the conjectures which have been described in this paper. 
It is unlikely that the method proposed here compares with the implicit technique 
of Douglas and Gallie [2] for computational efficiency, at least for one-dimensional 
problems. However, the present method seems well suited for computation in more 


dimensions. Furthermore, it is of interest to note that if, in obtaining Problem IIa, 
we were to replace q, by [¢?*’ — ¢77}]/Az, the resulting implicit difference equations 


would retain the feature of a sharply propagating interface and raise additional 
questions about the validity of our conjectures when extended to include implicit 
equations of this type. 
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Alternative Formulas for Osculatory and 
Hyperosculatory Inverse Interpolation 


By Herbert E. Salzer 


1. Introduction. In previous papers devoted to inverse osculatory interpolation, 
all of which were based upon the employment of the Lagrange-Hermite formulas 
for direct interpolation of f(z) (at equal intervals h for real functions or at points 
of a Cartesian grid of length h for complex functions) including both the osculatory 
(function with first derivative) and hyperosculatory (function with first and sec- 
ond derivative) types, the writer has given the Taylor series obtained by the in- 
version of the direct interpolation series about a suitable point z», in powers of 
r = [f(x) — f(ao)]/hf'(20) through the term in r” [1}-{3].* In fact, the same pro- 
cedure has been employed much earlier for inverse interpolation formulas based 
upon the ordinary non-osculatory Lagrange interpolation formula, which gives the 
inversion series in powers of a variable r that is proportional to f(x) — f(2») 
[4]+{6).* 

In this present note we give an alternative scheme for inverse osculatory or 
hyperosculatory interpolation in terms of the inverse function z(f), involving 
a(f) and dx(f)/df, or x(f), dx(f)/df and d’x(f)/df’ at separate points f; = f(z;). 
In other words, now we distinguish between the previously given point expansions 
in [1]}43], which may be characterized as special types of “inverse osculatory or 
hyperosculatory” formulas and the present “osculatory or hyperosculatory in- 
verse” formulas, which are analogous in structure to the osculatory or hyperoscula- 
tory direct interpolation formulas. 


2. Advantages of Alternative Scheme. 

A. The present scheme avoids any cumbersome explicit formulas by applying 
both a decomposability and uniqueness property of the Lagrange-Hermite inter- 
polation formula, which has been so effective in reducing the labor in direct inter- 
polation of high degree (described in detail in [3] and [7]), to those same kinds of 
osculatory formulas for the inverse functions. 

B. The formulas here, in terms of z(f) with either z’(f) or 2’(f) and 2x”(f) 
at f = f; , enable the user to go far beyond the 10th degree in accuracy with a frac- 
tion of the computational labor required for the power series formulas. Actual 
count of the number of operations for 10th degree accuracy by the older method 
and 11th degree accuracy (either 6-point osculatory or 4-point hyperosculatory ) 
using the present scheme, showed the latter to involve only around one-fourth of 
the number of operations required in the former. 

C. These alternative formulas are more truly interpolatory because of the ac- 
tual agreement with the inverse function and its derivatives at different points 
over the complete range, and the consequent closeness of the approximation near a 





Received November 20, 1959; in revised form, March 9, 1960. 

* This article is intended to be self-sufficient in the presentation of these alternative 
schemes for inverse interpolation. We shall avoid as much as possible the repetition of ma- 
terial in [1]-[8] to which the reader is referred for full details. 


257 











258 HERBERT E. SALZER 


number of points; whereas the convergence of the Taylor series expansion becomes 
rapidly poorer as | r | exceeds 1.* 

D. In this alternative scheme, since the f,’s are the fixed arguments and are 
unequally spaced, it hardly matters whether the corresponding z,’s are equally or 
unequally spaced. But for real variables, the previously given formulas in [1] and 
[3] would not be applicable when the z,’s are irregularly spaced. 

E. For complex osculatory or hyperosculatory inverse interpolation for z = z(f) 
when f(z) is an analytic function tabulated in the complex plane with f’(z), or 
with f’(z) and f”(z), there is no change in these present formulas other than the 
replacement of x by z. 


3. Osculatory and Hyperosculatory Inverse Interpolation Formulas. Before 
giving the alternative scheme for the osculatory and hyperosculatory cases, we 
should mention for the sake of completeness that even ordinary inverse interpola- 
tion has a corresponding alternative form. Instead of the formulas in [4]-(6], we 
may prefer the following concise rearrangement of Lagrange’s interpolation poly- 
nomial for the inverse function z(f): 


(1) olicss > axi/ dD ai, 


where 

(2’) a; = Ai/(f — fi), 
and 

(2”) A: = 1/ TI] (fs — fi). 


In the above (1), (2’), (2”), as well as in all subsequent formulas, x denotes either 
real or complex values, and in any >> or J the running index i, j or k has what- 
ever range is customary, say —[(m — 1)/2] to [n/2] for real interpolation or over 
any set of fixed grid points in complex interpolation (except for omissions indicated 
beneath the symbol). 

For all osculatory and hyperosculatory formulas we employ the first and second 
derivatives of the inverse function z(f), namely 2’(f) = 1/f’(x) and 2z”(f) = 
—f"(x)/[f'(z)}, at z = 2;, the same as at f(x) = f(z;), or more concisely, at 
f =f; . The notation z,’ and z,” is used for x’(f;) and x” (f;) respectively. It is not 
necessary to repeat here for the inverse function the development given in [3], 
[7] and [8] for concise expressions for the direct osculatory and hyperosculatory 
interpolation formulas, since those same ideas apply here. 

For the inverse function, in the osculatory and hyperosculatory formulas we 


* To see why, say in the real case, in a region where z(f) is one-valued and has derivatives 
of high enough order, consider the remainder term in the n-point osculatory or hyperosculatory 
Lagrange-Hermite formula for z(f), namely 


(1/(2n) 1) {TT Hn—aniay (F — $2) }2€d™2/(af)™) | sare 
or 
(1/(8n)!) (Ts 2fon-ayn F — fed) 2/(Af)™) | sary 


which can become very small for an f very close to f; , even though far from fy . 
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need the first and second derivatives of L;‘”(f), the n-point Lagrange interpolation 
coefficients in f defined by 


(3) Li (f) = i (f - i/T (f: — fi). 


Differentiation is with respect to f, after which we set f = f; .* These derivatives 
are conveniently ge as follows: 


(4’) = ae ODN bere = DUG: — fr, 
IF 
which is written more eu employing the notation 
(5) 1/(fi — fi) = aij, 
as 
(4) Lb" (fi) = Do ais. 
It 
Differentiating L;‘” (f) twice and setting f = f;, we obtain from (3), 
. fL.~ 
(6”) SPD =2 LDV - WG - We, 
7 /-_ a 


tke outside factor of 2 occurring because in every (j, k) combination there will be a 
(k, 7) combination, 7 ~ k. The double summation occurring in (6”) is avoided by 
employing the identity 2 }°>> aijaa = (>> a:;)* — } a7 so that 


(6’) Li" (fi) = (Qe as)’ — Deas, 
which from (4) is simply 
(6) LS" (Ff) = (Lief? - > a5 « 


For osculatory interpolation, following [7] p. 213, we define 
(7’) a; = A; 
which from (2”) and (5) may be expressed as 
(7) a; = iI a;}", 
and 
(8’) b; = —2A7L{""(f,), 
which from (2”), (4) and (5) is expressible as 
(8) b; = —2 {Tas}? 2 as,. 
From (7) and (8) we define 
(9) a: = a/(f — fi)’ + b/(f — fi), 





* For the occurrence of derivatives of Lagrangian coefficients in direct interpolation see 
(7] p. 213, for the osculatory case, and [3] p. 105, for the hyperosculatory case. 
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and 
(10) 8; = ai/(f — fi). 


Finally for the n points x; where we have both f(2;) and f’(z;), we find for z an 
approximation by the polynomial of the (2n — 1)-th degree in f = f(x), which 
is equal to x; at f = f; and whose derivative with respect to f is equal to 2) = 
z'(f:) at f = f;, according to 


(11) r~ > (aw; + Bw,’ />> Qi. 


For hyperosculatory interpolation, following [3] p. 105, and taking into account 
(2”), (4) and (5), we define 


(12) a; = {]] a;}*, 


IF#% 
€ 3 
Biff oul E eu. 
7#% 1-% 
and a quantity c; which in standard notation is given by 


(14’) ce; = As[—SL{""(f2 + 6fL{"(f,}"], 


(13) b; 


but when taking into account (2”), (4), (5) and (6), is expressible in present 
notation by 


(14) e: = ETT ss} ED as}? + $Y aw]. 


ivi ivi 

Next, using (12)—(14), we define 

(15) a: = ai/(f — fi)’ + b/(f — fi)’ + e/(F — fi), 
(16) B; = ai/(f — fi)” + bi/(f — fi), 

and 

(17) vi = a,/2(f — fi). 


Then for the n points x; where we are given f(x;), f’(z;) and f”(2;), we find for x 
an approximation by the polynomial of the (3n — 1)-th degree in f which, to- 
gether with its first and second derivative with respect to f, is equal to x;, 2x,’ 
and 2;” = —f"(2;)/[f'(x:)) at f = fi, according to 


(18) tT~ >> (ax; + Bie’ + vias”) / Do ai. 


In using (11) or (18), when there are many inverse interpolations with all 
the values of f being close to each other so that we have the same fixed points 
fi, the quantities a; and 6; in (7) and (8), or a; , b; and ¢; in (12)—(14) have to be 
computed just once, to be used repeatedly in (9)-(11) or (15)-(18). 


4. Application to Mathematical Table-Making. In a mathematical table where 
the inverse function will often be wanted, we might avoid the need for a separate 
table of the inverse function by supplementing the usual aids to direct interpolation 
(columns of differences or derivatives) with three, five, or even seven extra col- 
umns to facilitate osculatory or hyperosculatory inverse interpolation. For example, 
we might add three columns of just z;/ = 1/f’(z;) with a; and b; defined by (7) 








at 
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and (8) to aid just osculatory inverse interpolation by (9)—(11), or five columns 
of z/, 2,” = —f" (x) /f'(xoF, a;, 6; and c; defined by (12)-—(14) to help in just 
hyperosculatory inverse interpolation in (15)—(18), or even seven columns of z/’, 
z;”, one set of functions a; , b; defined by (7), (8) and another set of functions 
a;, b;, ¢;, defined by (12)—(14), giving the user a choice of either osculatory or 
hyperosculatory inverse interpolation. 

The use of these supplementary columns would not be restricted to tables of 
functions for just regularly spaced arguments z;. Thus we may tabulate these 
auxiliary quantities a;, b; and c; for tables having real arguments z; irregularly 
spaced or for tables having complex arguments in a Cartesian or polar grid. Even 
for osculatory or hyperosculatory direci interpoiation in a table whose arguments 
are irregularly spaced points z;, real or complex, if given f/ = f’(z;), or f/ and 
fi’ =f" (2), we may tabulate a; , b; by (7), (8), or a;, b; , e; by (12)-(14), and 
use (9)—(11) or (15)-(18), merely interchanging throughout in (5), (7)-(18) 
the variables z; with f; and xz with f. 

In the functions a; , b; and c; , the choice of the number n of fixed points should 
be, where feasible, sufficient to ensure full accuracy in the use of (11) or (18), 
which cannot be finer than the tabular uncertainty error of around ¢/f’(z;), the « 
being the error in the value of f(z;). 


Convair Astronautics 
San Diego, California 
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TECHNICAL NOTES AND SHORT PAPERS 


Permanents of Incidence Matrices 


By Paul J. Nikolai 


1. Introduction. This paper describes the evaluation of the permanents of cer- 
tain incidence matrices using the UNIVAC Scientific 1103A Computer. 

The permanent of an n by n matrix A = [a;;| with elements in a commutative 
ring is defined by the relation 


per (A) = dai, Arig *** Anig 
where the sum extends over all n! permutations of 1, 2, --- , n. The permanent is 
thus similar in definition to the determinant and suggests a theory for permanents 
analogous to that for determinants. No such theory is available, however, since 


permanents do not obey the analogue of the basic multiplicative law of deter- 
minants 


det (AB) = det (A) det (B). 
The analogue of the Laplace expansion is easily shown to remain valid for the 
permanent, but for computational purposes it is of little help. The calculations 
described in this paper were made practicable by a computational device given by 


H. J. Ryser. I shall state Ryser’s result, which has heretofore not been published, 
as a theorem. 


THEOREM 1. Denote by A, a matrix obtained from A by replacing r of the col- 


umn vectors of A by zero vectors. Denote by S(A,) the product of the row sums 
of A,. Then 


(1.1) per (A) = S(A) — 2) S(A1) + 2) S(A2) — «+» + (-1)"" BE S(Ana) 
1” 2” c” 
n—1 
where each sum extends over all the C,” ways of forming A, . 


Proof. The proof is based on the Principle of Inclusion and Exclusion or Sieve 
of Sylvester [2]. 


Consider an element of the form 
G@ = Ajj, Arje *** Anjn 


where the j; may assume any of the values 1, 2, --- , n. Let r denote the number of 
these integers not among ji, j2,---, jn . Suppose r > 0. Then a appears in (1.1) 


1 time in S(A), 
Cy" times in >> (Aj), 
C." times in >> S(A2), 


C,” times in > a S(A,) 
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or 


1-— C+ CC,’ —---+(-1)'C, =0 
times in all. 
On the other hand, if j; , jo, --- , j. is a permutation of 1, 2, --- ,n,r = O and 
a appears in (1.1) exactly once. 
One obvious advantage in using Theorem 1 is that the number of summands 


required to calculate per (A) is reduced from n! to 2”. Other features of this device 
are noted in Section 5. 


2. The v, k, \ Problem. Let it be required to arrange v elements into v sets such 
that every set contains exactly k distinct elements and such that every pair of sets 
has exactly \ elements in common, 0 < \ < k < v. This problem is referred to as 
the v, k, \ problem and the resulting arrangement is called a v, k, \ configuration or 
symmetric balanced incomplete block design. For a v, k, configuration list the ele- 
ments X, , X2,--- , X, in a row and the sets 7; , T: , --- , 7, ina column. Insert 1 
in row 7 and column j if X ; belongs to set T; and 0 otherwise. In this way is obtained 
a v by v matrix A of 0’s and 1’s called the incidence matriz of the v, k, \ configura- 
tion. It is not difficult to show that 


h = k(k — 1)/(v — 1) 
and that 


1 


| det (A) | = k(k — 2)? . 


Two v, k, \ configurations D,; and D, will be termed isomorphic if there is a one-to- 
one correspondence X, <> X2 = (X1)a between the elements {X;,} of D, and the 
elements {X-} of D. and a one-to-one correspondence 7; <> T,; = (71)8 between 
the sets {71} of D, and the sets {72} of De , such that if X, « T, , then (X;)a € (71) 8. 
An isomorphism of a design with itself is called a collineation. 

These combinatorial designs and their associated © .cidence matrices have been 
extensively studied. An excellent summary of the nzoblem together with an exten- 
sive bibliography can be found in [3]. 


3. Permanents of Incidence Matrices. A set R = {X,, X2,--- , X,} will be 
called a system of distinct representatives for the subsets T, , T2, --- , T, of the set 
D in case X;¢ T; and X; # X; fori # j,i = 1, 2,--- , n. The permanent of an 
incidence matrix possesses combinatorial significance in that it equals the number 
of systems of distinct representatives of the class of sets T,, T:,---, T,. As 
pointed out in Section 2, the determinant of the incidence matrix of a v, k, \ con- 
figuration is an elementary function of v, k, and \ alone. It seemed of interest to 
know whether or not the permanent possesses a similar property. Before trying 
calculations on UNIVAC Scientific, it seemed clear that per (A) could not be a 
simple function of v, k, and \ but ‘t remained an open question whether or not non- 
isomorphic designs having the same parameters v, k, and \ could possess unequal 
permanents. The answers to these questions would be of greater interest in the case 
of finite projective planes with n + 1 points per line which are », k, \ designs with 
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ven+n+i1,k n + 1, and A = 1. Unfortunately, the first case of non- 
isomorphic planes arises for n 9 with incidence matrices of order 91. Present 
theory does not permit calculation of the permanent of an incidence matrix of this 
size. There are known to be no instances of non-isomorphic v, k, \ configurations 
for v < 15. Different 15, 7, 3 designs, however, do exist, and the permanents of their 
incidence matrices were easily calculated on UNIVAC Scientific. 

Nandi [1] has constructed all 15, 7, 3 designs. There are five non-isomorphic 
examples. Study of these five designs revealed that exactly two which Nandi denotes 
by (vy’) and (a:a:’); possess a collineation of order 7 fixing one element (and one 
set) and permuting the remaining 14 elements (sets) in two cycles of length 7. 
For both designs let X; Ai, Az, --- , Ar; Bi, Be, --- , By denote the 15 elements. 
Each design has a set {Ai , Ao, --- , Az} fixed by the collineation (X) (A;A2 --- Az) 
(B, Bz --- Bz). (yy) and (a:a1’); can be displayed in the form 





{Ai, A2, As, As, As, Ao, Ay} 
{ X,Ai,A2, As, Bi, Be, By} 
{ X, Az, A;, As, Bo, Bs, Bs} 


{ X,A;,Ai, As, B;, Bi, Bs} 
{A,, Ar, As, Bs, Bs, Be , Bi} 
{A2, As, As, By, Be, B;, B;} 


{A;, Ai, As, Bs, Bs, Bs, Be 


and 


{Ay ’ Ao ’ A; ; Ag ’ As ’ Ags ’ A;} 


| X, 
{ X, 


{ X, 
{Ai, 
{Ao, 


{A7, 


Ai, Ac, As, Bi, Bs, Ba 
Ao ; Az ’ As ’ B, ’ By ; B;} 


A; ? A, ’ A; ’ B; 5] B:, B3} 


As ; Ag ’ B, ’ Bs ’ Bs 9 B;} 
A;, As, Bs, Be, Bz, By} 


Ai, Az, Bi, Ba, Bs , Be} 





respectively. Thus their corresponding incidence matrices appear as follows reflect- 
ing the collineation: 


011111110000000 011111110000000 
111010001101000 114010006 10110600 
101101000110100 101101000101100 
100110100011010 100110100010110 
100011010001101 100011010001011 
110001101000110 110001101000101 
101000110100011 101000111100010 
110100011010001 110100010110001 
0110106000010111 01104460001600111 
omer ee lt@@i001011 001101001010011 
000110101100101 0O00110101101001 
000011011110010 000011011110100 
010001100111001 010001100111010 
001000111011100 001000110011101 
010100010101110 010100011001110 


The permanent of each matrix is the sum of the equal minors belonging to the 
fixed row which represents the fixed set. This row and one column containing a 1 
in this row can be deleted reducing the order of the matrix by one. The permanent 
of the reduced matrix is computed and multiplied by the number of 1’s per row to 
yield the permanent of the incidence matrix. 


4. Results. The computational advantage offered by the collineations made 
the choice of (yy’) and (a:a’); for a first trial a natural one. The permanents of 
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the incidence matrices of these designs were found to be 24, 601, 472°and 24, 567, 
424 respectively. Thus nonisomorphic designs having the same parametets v, k, and 
\ may have unequal permanents, so that the permanent is not a function of », k, 
and } alone. 

The program for evaluation of these permanents turned out to be an efficient one 
requiring about four minutes of computer time for each matrix. The remaining 15, 
7, 3 designs were also run with the following results: 


Design Permanent 

(axyary")o 24 ,572 ,288 
(a2a") 24 ,567 ,424 
(6:61’) 24 ,582 ,016 


As computer time became available it was decided to calculate the permanents 
of the incidence matrices of all cyclic designs with prime v and v < 23. Cyclic designs 
with their simple structure and with v a prime might have yielded possible clues to 
a formula for the permanent in the cyclic case or to possible divisibility properties. 
Unfortunately no distinct cyclic designs with the same parameters arise for v < 31. 
Results of this survey are given as follows: 


uv k oN per (A) 

7 3 1 24 
7 4 2 144 
11 5 2 12,105 
11 6 3 75 .510 
13 t 1 3,852 
13 9 6 64 ,803 ,969 
19 9 4 142 ,408 ,674 ,153 
19 10 5 952 ,709 ,388 ,762 


In addition, the permanent of the incidence matrix of the projective plane of 
order 4, a cyclic 21, 5, 1 design, was computed and found to be 18, 534, 400. Here, 
as with two of the 15, 7, 3 designs, the use of a collineation reduced computation 
time by more than one half. This was a significant saving, better than four hours 
of computer time. 


5. Description of the Code for the UNIVAC Scientific Computer. Ones comple- 
ment binary arithmetic, two address logic, and an extensive array of logical instruc- 
tions together with equation (1.1) applied to 0,1 matrices contributed to a short, 
fast computer code for UNIVAC Scientific. The v rows of the square matrix A of 
0’s and 1’s were stored in the higher order v stages of v consecutive storage cells. 
An index r,0 S r S 2° — 1, counted the number of A,’s formed, and served as a 
logical multiplier in forming A, . If no row of A, were zero, S(A,) was calculated 
and added or subtracted from the accumulated sum according as the number of 
1’s in the binary representation of r was even or odd. A,_; was formed next and the 
calculation continued. The magnitude of the final sum was then per (A). 

All loops of the code were checked using the v by v matrix S of all 1’s. per (S) = 
v!. Machine accuracy was checked by running each calculation twice in every case. 
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On the Numerical Treatment of Heat Conduction 
Problems with Mixed Boundary Conditions 


By Arnold N. Lowan 


Abstract. The two-dimensional problem of heat conduction in a rectangle where 
the temperature is prescribed over a portion of the boundary while the temperature 
gradient is prescribed over the remainder of the boundary, may be treated nu- 
merically by replacing the differential equation of heat conduction and the equations 
expressing the given initial and boundary conditions by their difference analogs 
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Fig. 1.—Rectangular domain with ‘‘mixed’’ boundary conditions. 
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and solving the resulting system. It is shown that if the scheme is to be stable the 
intervals Az and Ay must be chosen so that kAt/(Az)* + kAt/(Ay)* < 3. 


Consider the two-dimensional problem of heat conduction in a rectangle, Figure 
1, when the temperature is prescribed over the thin line portion of the boundary, 
while the temperature gradient is prescribed over the heavy line portion of the 
boundary. This is a typical problem with “mixed” boundary conditions and should 
not be confused with the considerably simpler problem when the temperature is 
prescribed over certain complete sides of the rectangle, while the temperature 
gradient is prescribed over the remaining sides. As far as the writer is aware no 
analytical solution of the mixed boundary value problem above formulated (or 
of the analogous problem for the cylinder) is to be found in the literature. We must 
therefore (if interested in numerical answers) resort to the alternative of substituting 
for the differential equation of heat conduction and for the equations expressing 
the initial and boundary conditions their appropriate difference analogs, and solving 
the resulting system. 

The mathematical formulation of the problem is as follows: 


(1) oT (oT + ot) Osrsa0S5 y5bd,t>0 

(2) T(z, y, 0) = f(x, y) 

(3) T (0, y, t) = 0 0Sysb 

(4) T (z,0,t) =0 0Os2rtsGq 

(5) |2 T(z, y, o| = 0 as2zsa 

(6) |2 T(x, y, |, = 0 Osysb 
Ox =a 

(7) |2 T(x, y, »| = 0 cS2zsa 
oy => 

(8) T (z, 6, t) = 0 O0Osrs 


where for the sake of simplicity we have at first assumed that the prescribed tem- 


perature and temperature gradient are = 0. The difference analogs of the above 
equations are: 


(1*)) Trans = BT aaa +o Tin + (1 — 2a — 28) Tren 
+ a » |e +8 Thaticn h = i, 2, 3, viz M, k= l, 2, * N 


(2*) Tixo = f(hAz, kAy) 
(3*) Tern = 0 istah 
(4*) Tron = 0 1<h<o/Ar 


(5*) Tran = Th.o.n ¢,/ Az s h = a ‘Az 











268 ARNOLD N. LOWAN 


(6*) Tusten = Tukn lsk<N 
(7*) Tiwiin = Thowin @/Ar Sh S a/Ax 
(8*) Ti.wsin = 0 1 Sh S @/Ax 
where 

kAt kAt 


Taps. = T (hAz, kAy, nt) ; a= (Az)? ; B = (ay? 


At =a/(M+1) and Ay = b/(N +1). 


It will be convenient to consider the MN temperatures 7), ,,, withh = 1, 2,3, --- 
M and k = 1, 2, 3, --- N as the components of an M & N — dimensional vector to 
be denoted by T,,. It will also be convenient to replace the two subscripts h and k 
identifying the lattice point P, (i.e., the point with coordinate x = hAzx and y = 
kAy) by the single subscript p running from p = 1 to p = MN with the under- 
standing that for the M lattice points corresponding to k = 1, p runs from 1 to M, 
for the next set of M lattice points corresponding to k = 2 p runs from M + 1 to 
2M, --- etc. The system of MN equation (1*) may then be written in the matrix— 
vector form 


(9) Ta = AT. e 


As is well known, to prove the stability of (9), it suffices to prove that if S denotes 
the largest of the sums of absolute values of the elements of the rows of A then 
Ssl. 

Let 2 denote the set of lattice points closest to the boundary. It is clear that if the 
difference equation (1*) is applied to lattice points lying inside of ©, the resulting 
equation (which is in fact equation (1*) itself) has the five non-vanishing coeffi- 
cients 6, a, 1 — 2a —28, a and @. If we assume 1 — 2a —28 = Oora+8 S } 
it is clear that the sum of the absolute values of the coefficients is = 1. We shall 
show that if (1*) is applied to lattice points belonging to the set 2, the resulting 
equation may be characterized by the fact that the sum of absolute values of the 
coefficients is smaller than unity. Consider for instance the form taken by (1*) 
when applied to a point of 2 such that hAx < ¢, . Since for such a point Tho = 0 
the resulting equation has the four non-vanishing coefficients a, 1 — 2a —28, a and 
8. If again we assume 1 — 2a —26 > 0, it follows that the sum of the absolute 
values of the coefficients is = 1 — 8 < 1. In an entirely similar manner it is shown 
that if (1*) is applied to lattice points in 2 for which hAx = c, and the boundary 
condition (5*) is taken into account, the sum of the absolute values of the coeffi- 
cients of the resulting equation is 1 — a < 1. Similar conclusions may be drawn in 
the case of all lattice points in 2. Thus the quantity S previously defined is = 1. 
The stability of the difference scheme under consideration is thus proven, provided 
that the intervals Az and At are chosen so that 


kat kat 1 
et §@ 7-2 + 74; 


(Ax)? § (Ay)? ~ 2° 


In the above discussion we assumed that the prescribed temperature is = 0°C on 
the thin line portion of the boundary. If instead, nonvanishing temperatures are 
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prescribed on this portion, the criterion of stability is the same as before, since the 
error vector satisfies (1*) and evidently is = 0 on the portion of the boundary in 
question. 

We shall now discuss the modifications in the above analysis if the prescribed 


temperature gradient does not vanish. Let the above boundary conditions (5), (6) 
and (7) be replaced by 


oT 

(10) — = ¢ (z, t) a4 Srsay=0 
oy 
oT 

11 a = 

(11) = (y, t) z=a,0sysb 
oT 

(12) ay ~ & 8) @Szrtsay=b. 


The difference analogs of the last three equations are 


(10*) Tron = Tran + Ay-oi(hAr, nAt) o/At Sh S a/Az 
(11*) Tusien = Tun — At-d(kAy, nAt) lskswN 
(12*) Ti.wiin = Thwn — Ay-do(hAz,ndt) o/Ar Shs a/Az 


If the d:fference equation (1*) is applied to the lattice points for which the last 
three equations hold, we ultimately get 


(13) Trans = @Traan + (1 — 2a — 28)Taan + aT ryan 
+ BTron t+ Unan /Ae Sh S a/Az 
(14) Tan = BT wean + aT yin + (1 — 2a — 28)T wan 
+ BT urunnt Umen LSKEEN 
(15) Tawwntt = BT awn + Train + (1 — 2a — 28)Tiwin 
+ Trsiwin + Unwin 


where we have put 


Unan = BAy-di(hAz, nAt) c¢/Ar Sh S a/Az 
(16) Usain = —adzr-O(kAy, ndt) lskewN 
Unwin = — BAy-d2(hAz, nAt) Co/ Ax s h s a/ Az. 


The system of MN equations obtained by applying (1*) to the MN lattice 
points may be written in the form 


(17) Tau = AT, + U,, 


where U,, is a vector whose non-vanishing components are defined in (16) and whose 
remaining components are = 0. Since the error vector E, satisfies the difference 
equation (17) it is readily seen that 
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(18) Env = A"Ey + A” Up + A” UL +--- + U4. 


From (18) it follows at once that the criterion for the stability of (17) is identical 
with that for (9), namely that Az and Ay must be chosen so that a + 8 S }. 

It may be briefly mentioned that the above analysis may be extended to the 
more general case of the boundary conditions pT + 9q(d7T/dn) = F(t) where p and 
q take on prescribed values along the boundary. It may also be mentioned that the 
above analysis may be extended to problems with cylindrical and spherical sym- 
metry. 


Yeshiva University 
New York, New York: and 
AVCO Corporation 
Wilmington, Massachusetts 


High Precision Calculation of Arcsin x, Arccos z, 
and Arctan x 


By I. E. Perlin and J. R. Garrett 


1. Introduction. In this paper a polynomial approximation for Arctan x in the 
interval 0 S x S tan 2/24, accurate to twenty decimal places for fixed point rou- 
tines, and having an error of at most 2 in the twentieth significant figure for floating 
point routines is developed. By means of this polynomial Arctan z can be calculated 


for all real values of x. Arcsin x and Arccos x can be calculated by means of the 
identities: 


x : T 
Arctan Vi —_ Aresin 2 5 Arccos x 


2. Polynomial Approximation for Arctan x. A polynomial approximation for 
the Arctangent is obtained from the following Fourier series expansion, given by 
Kogbetliantz [1], [2] and Luke [3]. 


C) a ‘ 2i+1 
(2.1) Arctan (x tan 20) = 2 >> ( 1) (tan 6) Toi (x), 
i=0 2+41 
where 7';(x) are the Chebyshev polynomials, i.e., 7';(cos y) = cos (iy). The ex- 
pansion (2.1) is absolutely and uniformly convergent for |x| < land0 S @ < x/4. 


An approximating polynomial is obtained by truncating (2.1) after n terms. 
Thus, 





n-1 (__1)i 2i+1 
(2.2) P(x tan 20) = 2 > ( 1p Aten 6) Toi41(x). 





The truncation error is 


(2.3) | er | S tan 26 (tan 0)" | x |. 
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When z tan 2¢ is replaced by M and 7T2;4:(z) are expressed in terms of z, (2.2) 
becomes 


where 


k=0 k 


n—r—l 
B, = (1 — tan’9)”" >> (™ + *) (tan 0)™. 


With a choice of n = 9 and tan @ = tan 2/48, the following polynomial approxi- 
mation for Arctan z is obtained. 


(2.5) P(x) = ax + age’ +--+ + aye” 
where 
4% = 120 
a, = —0.33333 33333 33333 33160 7 
a; = 0.19999 99999 99998 24444 8 
a; = —0.14285 71428 56331 30652 9 
a@ = Q.11111 11109 07793 96739 3 


ay = —0.09090 90609 63367 76370 73 


a3 = 0.07692 04073 24915 40813 20 
a5 = —0.06652 48229 41310 82779 05 
&; = 0.05467 21009 39593 88069 41 


The truncation error is | er | < 6-10 | x |. 


3. Procedure for Calculation Arctan x. Subdivide the interval (0, ~) into seven 
intervals as follows:0 < u < tan 2/24, tan [(27 — 3)4/24] S u < tan [(2j — 1)x/24] 
for j = 2, 3, 4, 5, 6, and tan 117/24 < u < o. For | z| on the first interval use 
(2.5). For | x | on the (j + 1)st interval, (7 = 1, 2, 3, 4, 5), the formula employed 
is 


(3.1) Arctan | z| = 5 + Arctan t;, 
where 
~- wnt 
|a| — tan 12 
i; = ete oe 
1+ |2| tan 1D 


is used to obtain a value t; such that |¢;| < tan 2/24. Arctan ¢; is calculated by 
(2.5) and Arctan | x | from (3.1). When | z | is in the seventh interval 
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(3.2) Arctan |x| = . — Arctan Be: 
2 |x| 
and 
1 T 
eee - ee 
iz] = tan om 


The constants tan jr/24, (j = 1, 2, --- , 11) and x/2 are: 
tan 7/24 = 0.13165 24975 87395 85347 
tan 7/12 = 0.26794 91924 31122 70647 3 


bo 


tan 7/8 = 0.41421 35623 73095 04880 2 
tan 1/6 = 0.57735 02691 89625 76450 9 
tan 57/24 = 0.76732 69879 78960 34292 3 
tan 7/4 = 1.00000 00000 00000 00000 0 
tan 77/24 = 1.30322 53728 41205 75586 8 
tan 1/3 = 1.73205 08075 68877 29352 7 


tan 37/8 = 2.41421 35623 73095 04880 2 


tan 59/12 = 3.73205 08075 68877 29352 7 
tan 1lw/24 = 7.59575 41127 25150 44052 6 
a/2 = 1.57079 63267 94896 61923 1. 
4. Error Analysis. 


A. General. 

Errors arising from calculations by a computer may be classified into three 
categories according to Householder [4], namely: (1) truncation errors, (2) propa- 
gated errors, and 3) round-off errors. For the propagated error, if x and y are ap- 


proximated by x’ and y’, respectively, and the errors in each are denoted by ¢(2) 
and e(y), then: 


[e(a & y)| S |e(x)| + |e(y) |, 

















| ef 
jelav) | < |el2)| 4 ev) per 
ay x y 
| fx e(x (4 
\¢ (=) | “ ad ” 
i—— 1s — J v'y’ £0. 
= a e(y) 
: oe .- 7 
a . y 


For round-off error it is assumed that rounding is accomplished in the following 
manner. If \ digits are to b2 retained and the (A + 1)st digit is = 5, add one to 
the preceding digit; otherwise do not change the preceding digit. With this con- 
vention the round-off error in fixed point arithmetic is easily determined. For 
floating point arithmetic use is made of the following result. Let 





ee 
a- 
\p- 
x) 


yn- 
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x 


(a8 + af + --- + mB)8", - 0, 
(ysB* + yo + --- + w8)B", yn ~ 0, 


and z @ y represent addition, subtraction, multiplication or division. The round-off 
error in x @ y, «(x @ y), is given by 


|e(z @y)| S|z@y| (5) B*. 
Another result useful in floating point arithmetic is the following: If 


/ — 
~ 7 =| < af", (r = 1), then 2’ differs from z by at most a units in the rth sig- 


and y 





nificant digit. The preceding results are easily established. 
B. Errors in Arctan x. 


, B 
a) 0s 2 < tan ma" 

For fixed point arithmetic it shall be assumed that twenty-one decimals are 
used. The truncation error | er | < 8-10”. The error ¢, due to errors in the coeffi- 
cients in (2.5) is | ep| < 1.2-10°". The round-off error eg is | ez| < 6.35-10”. 
Hence, the total error is less than 8-10-”, and the calculated value of Arctan z is 
accurate to at least twenty decimal places. 

For floating point arithmetic using twenty-one significant digits, |er| < 
6-10 "x, | «>| < 10° ~z, and |ee| < 1.03-10°"z. Hence, | «(Arctan z)| < 
1.1-10° "2, and 


e(Arctan =) | 


< 12-10. 
Arctan x 


Thus the calculated value of Arctan zx differs from the true value by at most two 
units in the twentieth significant digit. 
b) tan <2 <tan DT. 

For fixed point arithmetic | ¢(t;) | < 7.8-10-”. The propagated error in Arctan 
t; due to this error in ¢; is | e(Arctan t;) | < 7.8-10-”. The total error in Arctan 1; 
is then | «(Arctan ¢;) | < 1.5-10-”. The error in Arctan z is then | e(Arctan x) | < 
2.5-10-", and the calculated value of Arctan z is accurate to twenty decimal places. 

For floating point arithmetic 

e(Arctan 2) | 





< 7.2-10™, 
Arctan x 
and the calculated value of Arctan x differs from the true value by at most eight 
units in the twentieth significant digit. 
llr 

c) tan=7 S2< 00. 

For fixed point arithmetic | «(Arctan x) | < 2.3-10™", and hence the calculated 
value is accurate to twenty decimal places. 
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For floating point arithmetic 


e(Arctan x) 


25-10” 
Arctan x <<a > 








and hence the calculated value differs from the true value by at most two units 
in the twentieth significant digit. 

C. Errors in Arcsin x and Arccos x. 

The error for floating point arithmetic using twenty-one significant digits will 
be given. Arcsin zx will be calculated by means of 


x 
Aresin x = Arctan rear t 


The quantity 1 — 2’ is calculated by means of 1 — 2 = (1 — x)(1 + x). Then 


e(Arctan x) 


. | 10°", 
; Aresin x ‘i 


and Arcsin z will be correct to within one unit in the nineteenth significant figure. 

The error for Arccos z is similar except that a round-off error due to subtraction 
is introduced. This error does not affect the conclusion that Arccos x will have been 
obtained correctly to within one unit in the nineteenth significant figure. 


5. Conclusions. From the standpoint of machine application the procedure given 
is economical and yields precise results. It uses only twenty stored constants; the 
calculation of Arctan x requires a maximum of only eleven multiplications and one 
division; the calculation of Arcsin x and Arccos x requires an additional multiplica- 
tion, division, and square root. 


Rich Electronic Computer Center 
Engineering Experiment Station 
Georgia Institute of Technology 
Atlanta, Georgia 


1. E. G. Koesetuiantz, “Computation of Arctan N for —» < N < using an electronic 
computer,’’ IBM Jn. Res. and Dev., v. 2, 1958, p. 43-53. 
2. E. G. Koasetiiantz, “Computation of Arcsin N for 0 < N < 1 using an electronic 
as 2s IBM Jn. Res. and Dev., v. 2, 1958, p. 218-222. 
. Y. L. Luxs, ‘On the computation of Log Z and Arctan Z,’’ MTAC, v. 11, 1957, p. 16-18. 
4. A. 8S. Housenouper, Principles of Numerical Analysis, McGraw Hill Book Company, 
New York, 1953. 


The Calculation of Toroidal Harmonics 
By A. Rotenberg 


1. Introduction. It is the purpose of this note to describe the mathematical tech- 
niques employed in a code [5] for the IBM 704 to calculate toroidal harmonics 
(associated Legendre functions of half-integral order). We use recurrence tech- 
niques similar to those used by Goldstein and Thaler [1] in calculating Bessel func- 
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tions. The functions of the first kind, P?_,(z), and the second kind, Q7_¢(z), both 
obey the recurrence relation 


(1) (vy — m)R,"(x) = (2v — 1)zRea(z) — (» + m — 1)R4(z) 


where we have written » = n — 3 and R,”(z) may be either P,”(z) or Q,"(z). Eq. 
(1) holds for any type of Legendre function. The toroidal harmonics are char- 
acterized by x = 1 and m and n integral. Since R7_4(z) = R™,4(xz), the computa- 
tions are carried out for n = 0 only. 

If two key values P,”(x) and P¥;:(x) are given, eq. (1) can be used to generate 
Pv+«(z) with good accuracy since P?,,(x) increases with, increasing k. However, 
the Q,”(x) decrease with increasing order and it is necessary to recur in the direction 
of decreasing order. The method of finding the key values is described in the next 
section. 


2. Method of Calculation. Many formulas exist for the calculation of P,”(z) 
[2], [3]. One which is well suited for machine computation is 


me. _ Tv +m+1) (2 — 1\""F(—»,m—»,m+1,(z—1)/(z + 1)) 
Se wy = wee (Ga) (e+ D> 


where F(a, b, c; z) is the hypergeometric function which can be computed from the 
series expansion [3]: 





. _) _ rte) $F TMa+n)r(b+n) , 
(3) F(a, b, 62) = FRG) 2, Tie+n)n! - * 


Since a and b are half-integral and c is an integer, the calculation is straightforward 
and the series converges very rapidly for all moderate values of xz. Since the first 
term in the series is unity and the series becomes monotonically decreasing after 
a finite number of terms, the computation continues until the next added term is 
less than some small value e. P™,; and P,;” are chosen as the key values and the re- 
currence relation is used to get functions of higher order. 

The calculation of the function of the second kind is somewhat more involved 
and a method suggested by Goldstein and Thaler [1] is used. Q7,,(z) is set to zero 
and Q,"(x) to a small value 6. Then the recursion relation is used to get the func- 
tions of lower order. It is necessary to begin computing with functions of fairly high 
order to get good accuracy. It is difficult to give an exact prescription, but for small 
values of x one should begin with the function of order 2» to get order » accurately. 
For values of z > 5, the functions increase very rapidly with decreasing order and 
only a few terms greater than v are needed. 

The values of Q,”(x) obtained by recursion are not properly normalized, and 
normalization can be effected by computing Q™,(z) say, from the formula [2]: 


Q,"(x) = VaI'(v + m + 1) ( P y" 1 


~— Qr#(y + 4)! z2—1 yt 











(4) 


2 ‘ 2 ss 
and normalizing the other computed functions to this value. For this normalization 
all values of Q,”(x) are positive. Many authors include a factor (—1)™ on the right 


Ast y—m+2 >t? 5) 
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side of eq. (4). The hypergeometric function appearing in eq. (4) can be computed 
as described above but convergence is slow for values of x close to unity, and in this 
case it may be desirable to compute the functions differently [4]. 

A subroutine [5] has been written for the IBM 704 to calculate the toroidal 
harmonics using the methods described here. For given m and z a table is obtained 
of the first twenty values of P?_4(r) and Q7_,(z) ie. n = 0, 1,--- , 19 as well 
as their derivatives er the formula [2]: 


(5) (2 - 1) 5 am ¢ Rz4(2) =(n — m + 3)Riu(x) — x(n + 4)RI4(2) 


Because the functions increase very rapidly with both m and z, it is convenient to 
make the restriction x < 40, m < 21. Where tabulated values exist [2], the code is 
found to give full agreement. In some other cases, the equation for the Wronskian 
[3] of the solutions was checked and found to be very accurately satisfied. The code 
computes correctly the toroidal harmonics to at least six significant figures. 

Atomic Energy Commission 

Computing and Applied Mathematics Center 


New York University 
New York, New York 

1. M. Goupste1n & R. M. Tuater, “Calculation of Bessel functions,’’ MTAC, v. 13, 1959, 
p. 102-108. 


2. NBS MartuematicaL TABLEs ey Tables of the Associated Legendre Function, 
Columbia University Press, New York, 


3. P. M. Morse & H. Fesnpacn, Methods of Theoretical Physics, McGraw-Hill Book 
Company, New York, 1953. 

4. J. P. Aurrray (private communication). 

5. A. Rorensera, NU ASLG, Share Distribution 711, Share Program Librarian, IBM, 590 
Madison Avenue, New York 22, New York. 


Transcendental Equation for the 
Schrodinger Equation 


By J. R. M. Radok 


The problem of the determination of the energy levels of single particles in 
cylindrical wells of different dimensions reduces to the transcendental equation 


h(i — Vie — &) + (ce) = 0 


for the Schrédinger equation, where ),"’, » are modified quotient Bessel functions 
for which a table has been published recently by Morio Onoe [1] and the variables 


2 2mUa ot _ ig, (E + Ua 
h2 ? h2 


involve the quantities 
m = mass, 
U = potential energy, 
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E = total energy, 
2ah 


Planck’s constant, 


a = radius of the well. 


Using this table, the first few roots have been obtained graphically and are 
recorded in Table 1 to three significant digits. For most practical purposes, these 


values should be satisfactory. If necessary, they can be improved by use of Newton’s 
method. 


Polytechnic Institute of Brooklyn 
Brooklyn, New York 


1. M. Onokg, Tables of Modified Quotients of Bessel Functions of the First Kind for Real and 
Imaginary Arguments, Columbia University Press, New York, 1958. 


A Note on Factors of ”' + 1 


By A. Gloden 


The factorizations enumerated in this note form a sequel to my published factor 
table [1] of integers n* + 1. They have been obtained by means of my table of 


solutions of the congruence z* + 1 = 0 (mod p) for primes lying between 
8-10° and 10° [2]. 


The following numbers are primes: 
n'+1 for n = 912, 914, 928, 930, 936, 952, 962, 966, 986, 992, 996. 
3(n' +1) for n 1071, 1087, 1101, 1119, 1123, 1125, 1135, 1163, 1173, 1183. 
e(n* +1) for n = 1562, 1726, 1732, 1834. 
Hi(n' +1) for n 


1818, 1848, 1982, 2006, 2012, 2064, 2088, 2094, 2228, 2340, 


2364. 
x(n’ +1) for n = 2346. 
2(n'+1) for n = 2262, 2302, 2544, 2682. 
ai(n' +1) for n = 2468. 
six(n* +1) for n = 2476. 


siz(n' +1) for n = 2808. 








sig(n' +1) for n = 1709, 1715, 1759, 1787, 1827, 1845, 1855, 1879, 1895, 1963, 
2015, 2021, 2031, 2093, 2185, 2229, 2259, 2287, 2303, 
2327, 2331. 

saq(n' +1) for n = 2211, 2299, 2651, 2761, 2791, 2815. 


siag(n' +1) for n = 2533, 2577, 2691, 2723, 2857. 
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st5(n' +1) for n = 2747, 2771, 2885. 
s4q(n' +1) for n = 2669, 2683, 2749. 


New factorizations are as follows: 


938° + 1 = 809273-956569 
1060 + 1 = 847577-1489513 
1348 + 1 = 940169-3511993 
1512° + 1 = 926617-5640361 
1874 + 1 = 914561-13485457 
2100‘ + 1 = 17-873553-1309601 
2838' + 1 = 868841-74663657 
2908‘ + 1 = 41-940369- 1854793 


4(1155° + 1) = 830233-1071761 
4(1191‘ + 1) = 935353-1075577 
4(1509° + 1) = 872369-2971849 
4(2635* + 1) = 857569-28107577 
4(2765° + 1) = 908353-32173321 
4(2977* + 1) = 17-809041 -2855393 
The following factorization was omitted from my original table [1]: 
4(2055* + 1) = 17-572233-916633. 


The least integers still incompletely factored correspond ton = 1038 and 1229, 
for even and odd values of n, respectively. 


11 rue Jean Jaurés 
iio 


A. GiopeEn, ‘‘Table de factorisation des nombres n‘ + 1 dans |’intervalle 1000 < n < 
3000, m Tastieut Grand-Dueal de Luxembourg, Archives, Tome XVI, Luxembourg, 1946, p. 71- 
88. 


2. A. GLopEn, Table des Solutions de la Congruence x* + 1 = 0 (mod p) pour 800,000 < p 
< 1,000,000, published by the author, rue Jean Jaurés, 11, Luxembourg, 1959. 


A Note on the Solution of Quartic Equations 


By Herbert E. Salzer 
For any quartic equation with real coefficients, 
(1) X* + AX* + BX’? +CX+D=0, 


the following condensation of the customary algebraic solution is recommended as 
quickest and easiest for the computer to follow (no mental effort required). It 
works in every exceptional case. 
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Denote the four roots of (1), by Xi, X2, X;, and X,. With the aid of [1], solve 


the “resolvent cubic equation” az* + bz” + cx + d = 0 for the real root x, only, 
where 


(2) a=1, b=-B, c=AC-—4D, and d= D(4B-— A’) -—C. 
Find 


(3) m= ++V/1tA? —-B+n, n= Ax, — X 
4m 
If m = 0, take n = +/ix,;2 — D and proceed according to the following Case I or 
Case II, depending upon whether m is real or imaginary. 
Case I: If m is real, let ($A — 1, — B) = a, 4n — Am = B, Va + B =7, 
Va — B = 4, and finally 








ou an aie 
oN I ee — Sol 2 oF 
(41) , yi 
dee _ am 
X3 = 44 +” y ? and Xs = 2A ots $ 





Case II: If m is imaginary, say m = im’, then n is also imaginary, say n = in’. 
Let 


(3A? — 2 — B) =a, 4n’ — Am’ = 8B, +Vo? + # =p, 


and finally 


x, = Att i(m' + 8) 
2 ? 





X; = X,, the complex conjugate of X,, 
(411) 


_ _ fA —y+ilm'—3 
xX, = — i 








and X, = X;, the complex conjugate of X;. 


If y = 0, we must have a = —a’, a’ = 0, and formula (4II) still holds provided 
that in it we replace 6 by + Va’. 

As an example consider the quartic equation X* + X*° + X°+ X¥ +1 =0, 
where A = B = C = D = 1, so that from (2) the resolvent cubic equation is 
x —x —3c+2=0. From [1] we find 2, = 0.61803 400. From (3), m = 
++/-=0.13196 600 = +0.36327 125i, so that m’ = +0.36327 125. Then n = 
128196 OD nn = +0.95105 6557, so that n’ = +0.95105 655. Proceeding according 

1.45308 5002 , ue : 

to Case IT, a = —1.11803 400, B = 3.44095 495, p = 3.61803 41, y = 1.11803 40 
and 6 = 1.53884 18. Then from (4I1) we obtain X, = 0.30901 70 + 0.95105 65:, 
X, = X, = 0.30901 70 — 0.95105 657, X; = —0.80901 70 — 0.58778 53% and 
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X, = X; = —0.80901 70 + 0.58778 53i. These roots may be verified as correct, 
since they are known to be the complex fifth roots of unity, namely X, = cos 72° + 
i sin 72°, X: = cos 288° + i sin 288°, X; = cos 216° + i sin 216°, and X, = 
cos 144° + 7 sin 144°. 

Convair Astronautics 


San Diego, California 


1. H. E. Savzer, C. H. Ricnarps & I. Arsuam, Table for the Solution of Cubic Equations, 
McGraw-Hill, New York, 1958. 


A Conjugate Factor Method for the Solution 
of a Cubic 


By D. A. Magula 


1. Introduction. This paper gives a simple method for computing the real roots 
of the reduced cubic equation with real coefficients, 


(1) z+Azr+Be=0, 
having roots a, b, c. We assume a to be real, since every cubic equation has at least 
one real root. 


The method consists in factoring B, and setting one factor equal to ++/m, 
the other n. For all pairs m, n such that m + n = —A, ++/misa root. If no such 
pair exists, a method of interpolation is shown. 


2. Proof of Method. The reduced cubic equation (1) can be transformed, by 
using the relations between the roots and coefficients, into a complete cubic, 





(2) p + 6Ap’ + 9A*p + 4A* + 27B’ = 0, 
where 

(3) p = (—3a" — 4A). 
Equation (2) can be written in the form: 

(4) (p + A)*(—p — 4A) = 27B° 
or 

(5) (p as #2 o/ > = oe we. 
Let 

(6) a0 —_— 7 ¢ and n= Po 4 
(7) n/m = +B 

and 

(8) m+n= —A. 





Received September 21, 1959; in revised form, December 22, 1959. 











282 D. A. MAGULA 
From (7) and (8) 
Vm(—A —m) = +B 
or squaring both sides: 
(9) m* + 2Am’ + A’m — B’ = 0. 


But the equation (9) is identical with the equation (1) whose roots are squared. 
Therefore, 


(10) m=a. 


Thus one root of the given cubic is ++/m. The other two roots can be found by 
ealving the quadratic equation resulting from the depression of the given cubic: 


(11) e+art+ (A + a) = 0. 


3. Steps in Computation. The calculation can be carried out in the following 
eiraple operations. A model example is given: «° — 7x — 6 = 0. 

A. Make a table with pairs of integral factors of B putting one in the first 
column, +~+/m, and the other in the second, n. In each case where m + n = —A, 
the corresponding entry in the first column is a root. 




















+ Vm n m m+n 
1 —o 1 —5 
2 —3 4 1 
® -2 9 | 7 
6 —] 36 35 
-@ 6 1 7 
-@ 3 4 | 7 





In this example the equation has three integer roots which appear circled in the 
first column. 

B. If the equation has one integer root, it can be found by the above method. 
Then the other two can be found by solving equation (11). 

C. If the equation has no integer roots, we may use this method to interpolate. 
For example, suppose we have the equation, z* — 3 = 0, in which A = 0. We find 
from the corresponding table that the root is between 1 and 3. We can then make 
additional entries in the table by taking values within this interval. 




















+VJVm n m m+n 
1 —3 1 —2 

3 —i1 9 +8 

1.4 —2.14 1.96 — .18 
1.5 —2 2.25 + .25 
1.44 —2.08 2.07 — .01 
1.45 —2.07 2.10 + .03 
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Thus the real root of this equation is between 1.44 and 1.45, and*this process 
can be repeated until the desired degree of accuracy is attained. From (B) it 
follows that the presence of imaginary roots in no way affects the solution of the 
cubic by this method. 


817 West End Avenue 
New York, 25, New York 
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42(F].—A. Giopen, Table des Solutions de la Congruence x* + 1 = 0 (mod p) pour 
800 000 S p S 1000000, published by the author, rue Jean Jaurés, 11, 
Luxembourg, 1959, 22p., 30 cm., mimeographed. Price 150 Belgian francs. 


This volume represents the culmination of independent efforts of table-makers 
such as Cunningham, Hoppenot, Delfeld, and the author, extending over a period 
of several decades. The foreword contains an extensive list of references to earlier 
publications of this type, which combined with the tables under review give the 
two least positive solutions of the congruence x* + 1 = 0 (mod p) for all admissible 
primes less than one million. [See RMT 109, MT AC, v. 11, 1957, p. 274 for a similar 
list of references.] 

Professor Gloden has used such congruence tables in the construction of manu- 
script factor tables of integers N* + 1, which now extend to N = 40000, with 
certain omissions. The bibliography in the present set of tables also contains refer- 
ences to this work. Numerous references to these factor tables are also listed in 
RMT 2, MT AC, v. 12, 1958, p. 63. 


J. W. W. 


43(G, X].—F. R. Gantmacuer, Applications of the Theory of Matrices, translated 
by J. L. Brenner, Interscience Pub., New York, 1959, ix + 317 p., 24 em. 
Price $9.00. 


This is a remarkable book containing material, not easily available elsewhere, 
related to the analysis of matrices as opposed to the algebra of matrices. In this I 
use the word analysis to mean broadly that part of mathematics largely dependent 
upon inequalities (and limits) as opposed to algebra, which depends largely on 
equalities. 

In particular, the material in this book is directed largely toward studies of 
stability of solution of linear differential equations (in Chapters IV and V) and of 
matrices with nonnegative elements (in Chapter ITI). 

Several aspects of the book will be useful to numerical analysts. These include 
some parts of the chapter on matrices with nonnegative elements, the implications 
of the chapters on stability of solutions of differential equations to the stability of 
numerical methods of solving differential equations, and (as the author points out) 
a numerically feasible method of finding the roots of polynomials. 

Many topics included in this text are not easily available elsewhere. For example, 
product integration is expounded; this is an amusing version of Euler’s method 
applied to the solution of linear first-order homogeneous differential equations. 
Most of the other exposition is unique in one way or another, and on the whole the 
book is a valuable contribution to literature. 

This is a translation, augmented to some extent by bibliographic and other 
notes, of the second part of a successful Russian book. It is interesting to note 
that another publisher has announced the impending publication of translations of 
both parts. 

The printing is good, and the reviewer noticed no serious errors. There are four 
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short appendices devoted to standard elementary theorems, a bibliography with 
seventy-two entries, and a useful index. The chapter headings follow: 

Chapter I. Complex Symmetric, Antisymmetric and Orthogonal Matrices 

Chapter II. Singular Bundles of Matrices 

Chapter III. Matrices with Nonnegative Elements 

Chapter IV. Applications of the Theory of Matrices to the Study of Systems of 

Linear Differential Equations 
Chapter V. The Routh-Hurwitz Problem and Related Questions. 
6.2 


44/K).—Joseru Berkson, “Tables for use in estimating the normal distribution by 
normit analysis,” Biometrika, v. 44, 1957, p. 411-435. 


In a quantal response assay a number of independent tests are made at each of a 
number of dose levels, and the result of each test is graded as “‘success” or “failure’’. 


If the probability of “success” at dose metameter value z is assumed to follow the 
“normal” law 


P( n (z—#) /¢ re 
x) = 1/V2e € dt 


the method of “normit analysis” is proposed by Berkson as a replacement for the 
familiar (iterative) method of “probit analysis” for estimating the parameters 
u and o. 

Suppose that the dose metameter values used are x , --- , 2 , that n; tests are 
made at level x; , and that r; of these tests result in success. Let 


1/(2n;) if T= 0 


pi = 4ri/niffO <r< nn; 
1 — 1/(2n,) if r; = nj, 
X; = X(pi), where X(p) is defined by the relation 
X(p) ’ 
p = (1/V/ 2x) / e”” du, 
Zi = (1/V/2n) *"", 
a = Z?/pdl Fis pi). 


The method consists of a weighted regression analysis, which is facilitated by tables 
which give for each p; the corresponding values of w; and w;,X; . 

In table 2 w; and w;X; are given to 6D for p = 0.001(0.001)0.500. (For p > 3, 
w(p) = w(l — p) and wX(p) = —wX(1 — p).) For moderate n; interpolation 
may be avoided by use of Table 1, which gives w; and w;X;, also to 6D, for all 
combinations of r; and n; for which 1 < n; S 50 and 0 S r; S n/2. (Forr > n/2, 
w(r,n) = w(n — r,n) and wX(r,n) = —wX(n — 1, n).) It is stated that the 
entries in both tables are correct to within +1 in the final digit. 


Paut MEIER 
University of Chicago 
Chicago, Illinois 
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45(K).—I. D. J. Bross & E. L. Kasten, “Rapid analysis of 2 x 2 tables’, Amer. 
Stat. Assn., Jn., v. 52, 1957, p. 18-28. 


Conventional statistical analysis of 2 x 2 tables such as 











Sample A A 
1 a b N, = NP 
2 | c d Nz = NQ 
| 1 N 





involves use of triple-entry tables for critical values of a. These tables are entered 
with N, Ni, and 7, or some equivalent combination of three numbers. The body 
of the table then usually gives critical values for the observation a. The authors 
remark that the statistical test is relatively insensitive to variation in N and pro- 
pose to reduce the complexity of the tabular entry to double entry by ignoring N 
and using only the parameters 7 and P. Charts I to IV inclusive present lower tail 
critical values for a at 5%, 2.5%, 1% and .5% levels of probability for .1 <P < .9 
and 5 = T < 49. Interchange of P and Q produces lower tail critical values for c 
(and by subtraction) upper tail critical values for a. 

The authors claim that the approximation is good, provided P and Q are both 
at least .1 and T is not larger than .2N. 


Leo Katz 
Michigan State University 
East Lansing, Michigan 


46(K].—F. E. Cuark, “Truncation to meet requirements on means,’’ Amer. Stat. 
Assn., Jn., v. 52, 1957, p. 527-536. 


The problem under consideration is that of truncating a given distribution so 
that the resulting population will meet specified sampling requirements. This 
problem arises when one wishes to screen the output of some production process in 
order to reduce the risk (probability) of having lots rejected on the basis of a re- 
quirement that only those lots will be accepted for which the mean X of a random 
sample of n items shall, for example, exceed or be less than some value, say UAL 
(upper average level) or LAL (lower average level). 

Methods are given for determining a single point of truncation A such that the 
mean X, of a random sample from a normal population (u, «) screened or truncated 
at A will meet a specification requirement X, = LAL or X, < UAL with risk of 
rejection r. 

Methods are also given for determining double points of truncation A and B 
such that a normal population (4, o) truncated at X = A and at X = B will meet 
the requirement LAL < X42 S UAL with risk r. 

As aids in carrying out the computations involved in the above methods, a table 
is included which lists values of the mean p,, and the standard deviation o,, of 
the standard normal population (0, 1) truncated at a and at b (a < b). Entries are 
given to 4D fora = —3.00(.25).50 and for b = 3.00(—.25)0. A chart is included 
which contains curves of constant ya, and o, for fixed degrees of truncation p, 
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where p is the proportion of the complete population which is eliminated’ by trunca- 
tion. In this chart, a extends from —3.0 to 0.5, b extends from —0.5 to +3, p = .05, 
-10(.10)1.0, was = —1.0(.1)1.0, a4, = 0(.1).9. A second chart contains a set of 
five curves for selected values of n and r to be used in determining a and p as a 
function of h, wnere h = («4 — LAL)/c. Values of a extend from —1.4 to 0.4, h 
extends from —1.0 to 0.4, and p from .10 to .65. 


A. C. Conen, Jr. 
University of Georgia 
Athens, Georgia 


47(K].—W. H. CLatwortuy, Contributions on Partially Balanced Incomplete Block 
Designs with Two Associate Classes, NBS Applied Mathematics Series, No. 47, 
U. 8. Government Printing Office, Washington 25, D. C., 1956, iv + 70 p., 
26 em. Price $.45. 


This publication contains six papers dealing with various aspects (enumeration, 
dualization, and tabulation) of partially balanced incomplete block designs with 
two associate classes, and with the construction of some new group divisible designs, 
triangular incomplete block designs, and Latin square type designs with two con- 
straints. Approximately 75 new designs not contained in the monograph of Bose, 
Clatworthy, and Shrikhande [1] are given in the present paper. A number of theo- 
rems are proved in the six papers. Two of the theorems give bounds on the parame- 
ters v, pi , and py in terms of the parameters r, k,n; , nz , 1 , and ds of the partially 
balanced incomplete block design with two associate classes. The two theorems 
on the duals of partially balanced designs are useful in identifying certain partially 
balanced incomplete block designs. 


W. T. Feperer 


Cornell University 
Ithaca, New York 


1. R. C. Boss, W. H. Ciratwortny & S. 8. SurikHanpe, Tables of Partially Balanced 
Designs with Two Associate C lasses, North Carolina Agricultural Experiment Station Technical 
Bulletin No. 107, 1954. 


48/K|.—W. J. Drxon, “Estimates of the mean and standard deviation of a normal 
population,” Ann. Math. Stat., v. 28, 1957, p. 806-809. 


Four estimates of the mean in samples of N from a normal population are com- 
pared as to variance and efficiency. These are (a) median, (b) mid-range, (c) mean of the 
best two, (d) Xy.wc = > oN [X./(N — 2)]. The sample values are denoted X, < 
X. S --- S Xy. The results for the median and mid-range are given primarily for 
comparison purposes, since results are well known. The mean of the best two is 
reported as the small sample equivalent of the mean of the 27th and 73rd per- 
centiles. 

The variance and efficiency are given to 3S for N = 2(1)20. The estimate (d) 
is compared to the best linear systematic statistics (BLSS) as developed in [1] 
and [2]. It is noted that the ratio Var (BLSS)/Var (X),,v;) is never less than 0.990. 

Two estimates of the standard deviation are given in Table II. One, the range, 
is well known. The quantity k which satisfies E(kW) = o is tabulated to 3D for N = 
2(1)20. Denote the subranges Xy_ii, — X; by W,) and Wa) = W. The estimate 
S' = k'(>>W.), where the summation is over the subset of all W,,) which gives 
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minimum variance, is the other estimate. Table II compares variances and efficiences 
of these two estimates to 3S for N = 2(1)20. Also a column gives the ratio of the 
variance of (BLSS) as given in [2] to the variance of S’ to 3D for N = 2(1)20. 


J. R. VaTNsDAL 
State College of Washington 
Pullman, Washington 


1. A. K. Gupta, “Estimation of the mean and standard deviation of a normal population 
for a censored sample”’, Biometrika, v. 39, 1952, p. 260-273. 

2. A. E. Sarwan & B. G. GREENBERG, “Estimation of location and scale parameters by 
order statistics from singly and doubly censored samples. Part I’’, Ann. Math. Stat. v. 27, 
1956, p. 427-451. 





John Wiley & Sons, Inc., New York, 1959, xi + 224 p., 29 cm. Price $8.00. 


The first feature one notices about the second edition of the Dodge-Romig 
tables, as compared to the first edition, is its size. Measuring 8}” x 11”, with a total 
of 224 pages, it stands in bold contrast to the pocket-sized 53” x 84”, 106-page first 
edition,—a four-fold increase, in a sense indicative of the growth of statistical 
quality control since the first edition was published in 1944. 

There are 60 pages of text covering an introduction, and four chapters which 
describe the principles, procedure for application, and the mathematics of the 
sampling plans. The titles of these chapters are, respectively: A Method of Sampling 
Inspection; Single Sampling and Double Sampling Inspection Tables; Using Double 
Sampling Inspection in a Manufacturing Piant; and Operating Characteristics of 
Sampling Plans. The remaining 158 pages include a table of contents, seven ap- 
pendices, and an index. The last four appendices give the same four sets of tables 
as appear in the first edition; these are respectively entitled: Single Sampling 
Tables for Stated Values of Lot Tolerance Per Cent Defective (LTPD) with Con- 
sumer’s Risk of 0.10; Double Sampling Tables for Stated Values of Lot Tolerance 
Per Cent Defective (LTPD) with Consumer’s Risk of 0.10; Single Sampling Tables 
for Stated Values of Average Outgoing Quality Limit (AOQL); Double Sampling 
Tables for Stated Values of Average Outgoing Quality Limit (AOQL). The first 
three appendices are devoted to 120 pages of operating characteristic curves; these 
are, respectively, OC Curves for all Single Sampling Plans in Appendix 6; OC Curves 
for all Double Sampling Plans in Appendix 7; and OC Curves for Single Sampling 
Plans with c = 0, 1, 2, 3, and m S 590 (based on binomial probabilities). 

As can be seen from the above list of contents, the largest part of the increase 
in size is due to the inclusion of three sets of operating characteristic curves,—two 
sets for the AOQL plans, and one set for a separate inventory of single sampling 
plans. This separate inventory has a wide enough range in sample size and ac- 
ceptance numbers to provide a useful independent reference of OC curves for those 
who prefer to derive their own single sampling plans (charts to derive such plans 
have been retained from the first edition). Although separate sets of OC curves are 
not given for the LTPD plans, the authors point out how these may be obtained 
or estimated from the OC curves which are given. 

It is heartening to see this inclusion of OC curves and the authors’ comment 
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that it “has been urged over the years by a number of engineers”, as well as the 
inclusion of a completely new chapter devoted to a discussion of the operating 
characteristic curves. This overt endorsement by this distinguished team in the 
field of quality control of the use of the OC curve to evaluate an acceptance sampling 
plan should impress upon quality control practitioners the importance of describing 
the assurance provided by a decision-making procedure in probability terms. Per- 
haps it will also emphasize the fact that while one may prefer to attach the label of 
LTPD plan, AOQL pian, or AQL plan to an acceptance sampling plan, depending 
on first considerations in deriving or classifying a plan, these sampling plans are 
one and the same as far as assurance in decision making is concerned, if they have 
the same operating characteristic curve. 

The clear distinction made by the authors between OC curves giving the prob- 
ability of lot acceptance based on lot quality as distinguished from process quality 
is most welcome, as this distinction is seldom clearly made. That there is a difference 
is often not realized; sometimes it is misunderstood and the importance of the 
difference exaggerated; at best it is ignored, since most often, but not always, the 
difference in OC curves is slight, as the authors point out. It is, however, unfortu- 
nate, in this reviewer’s opinion, that the authors chose to attach the special labels 
of Type A and Type B to the corresponding OC curves, rather than simply to 
identify them as finite lot quality and infinite lot or process quality OC curves. The 
add‘tional labels are not essential, and can do nothing more than add mystery to 
an already confused situation to the many who will not look beneath the labels. 
Unfortunately, also, an inaccuracy has crept into a statement with regard to these 
OC curves. On page 59 starting at the bottom of the first column the authors state, 
“When the sample sizes are a larger percentage of the lot size, the Type A OC curve 
will fall somewhat below the Type B curve shown on the chart, as can be seen in 
Fig. 4-1 where the Type A OC Curve for N = © is identically the Type B OC 
curve.”’ That the statement is inaccurate can indeed be seen from a careful examina- 
tion of Fig. 4-1, since the finite lot (Type A) and infinite lot (Type B) OC curves 
intersect and cross. It would be better to remember that for the same sampling plan 
the OC curve for a finite lot is always more discriminating than the OC curve for an 
infinite lot. 

Little need be said about the tables of sampling plans, which are well known as 
a result of the first edition. Derived on the principle of minimizing the total amount 
of inspection, sampling as well as screening of rejected lots, they are particularly 
suited to producers who are responsible for both the production and sampling inspec- 
tion of their finished product. The sampling plans may also be used by a purchaser 
for acceptance inspection, but the choice of plan for this purpose should be based 
principally on the properties of the operating characteristic curve. 

The format and typesetting, including tables and graphs, are considerably 
improved over the first edition. The division of each page of text into two columns 
also improves on readability. A minor typographical error which occurs in the 
second edition, but not in the first, appears on page 33, equation (2-la), where 
Cy should be CX. 

The book can be highly recommended to those with modest or little mathemati- 
cal background. The improvements in the second edition are sufficient to warrant 
its own place, along with other worthy texts, on the bookshelf of students and 
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practitioners of quality control who are interested in a comprehensive account of 
sampling inspection as well as in the procedures and tables for its application. 


Harry M. Rosensiatr 
Federal Aviation Agency 
Washington, District of Columbia 


60(K].—C. W. Dunnett & R. A. Lamm, “Some tables of the multivariate normal 
probability integral with correlation coefficients 4,’’ Lederle Laboratories, Pearl 
River, New York. Deposited in UMT File. 


The probability integral of the multivariate normal distribution in n dimensions, 








having all correlation coefficients equal to p (where necessarily — +5 <p<1l), 
is given by 

bn ua r (¢ 1 y" (1 + (n — A)” ma [ 1 + (n al 2)p_ % 

a ee Go aeee PLO = pl + @ — ipl 


\ coc ie 
{= “ “ite “Bp 2y Pe Bip Jan —— 


This function, which we shall denote by F,,,0(2 , --- , kn), has been tabulated for 
p = 4 and x, = --- = 2, by Teichroew [1]. In the present paper, we present a 
table for the case p = 34 and 2, = --- = z,. The need for this table arose in con- 
nection with a multiple-decision problem considered by one of the authors [2]. 

In computing the table, use was made of the fact that, for p = 0, Fn.p(m,--- , 
Zn) belongs to a class of multivariate normal probability integrals which can be 
written as single integrals (see Dunnett and Sobel [3]), a fact which greatly facili- 
tates their numerical computation. In this case, we have 


+o 2 -, 
ae p(x, is Saat In) = / II [ (EF) | 5 dy 
—o i=] —@ 





where 


fly) = es “’* and F(y) = [ f(y) dy. 


The attached table was computed by replacing the right-hand side of this identity 
by the series based on the roots of Hermite polynomials described by Salzer et al. 
[4]. Those tabular values marked with an asterisk have been checked by comparison 
with the values obtained by applying Simpson’s rule. The values checked were 
found to be systematically less than the Simpson’s rule values by an amount which 
varied between .0000000 and .0000013, depending on n. This indicates that the error 
in the tabular values may be no more than 1 or 2 units in 6th decimal place, but 
further checks are required in order to substantiate this. 

The table gives F,, 1/3(x, --- , x) to six decimal places, with x varying from 0 to 
7.0/+/3 in steps of 0.1/+/3 for n = 1 (1) 10, and from 1.5/+/3 to 2.1/+/3 in 
steps of 0.01/+/3 for n = 1 (1) 10, 13, 18. 
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1. D. Te1curoew, ‘Probabilities associated with order statistics in samplés from two 


normal populations with equal variance,’’ Chemical Corps Engineering Agency, Army Chemi- 
cal wk Maryland, 1955 


. DUNNETT, “On waeriee by 1000 largest of k normal population means,’’ (to be pub- 
lished in J n., Roy. Stat. Soc, Series . 


_ C. W. Dunner & M. ‘Soe, “A proximations to the probability integral and certain 


Ren ae points of a multivariate analogue of Student’s t-distribution,’’ Biometrika, v. 42, 
rere p. 258. 


H. E. Sauzer, R. Zucker & R. Capuano, ‘“Table of the zeros and weight factors of the 
dest | twenty Hermite polynomials,”’ Jn. Res., Nat. Bur. Standards, v. 48, 1952, p. 111. 


§1(K].—E. C. Frevier, H. O. Hartitey & E. 8. Pearson, “Tests for rank correla- 
tion coefficients. I,”” Biometrika, v. 44, 1957, p. 470-481. 


This paper is concerned with sampling determination of the approximate distri- 
bution for zs = tanh” ’rs and zx = tanh” 'rx , where rs is Spearman’s rank correla- 
tion coefficient and rx is Kendall’s rank correlation coefficient, for the case of sample 
of size n from a bivariate normal distribution. It is concluded that zs and zx are 
approximately normally distributed if n is not too small, with var (zs) = 1.060/ 
(n — 3) and var (rg) = 0.437/(n — 4). Eight tables are presented. Table 1 con- 
tains 4D values of three versions of var (rs) for p = 0.1(0.1)0.9 and n = 10, 30, 50; 
one version is Kendall’s approximate formula (adjusted), another is the observed 
value, and the third is a smoothed form of the observed value. Table 2 contains 
3D values of var (rs)/{1 — (Ers)*| and 4D values of var (rx)/{1 — (Erx)’|, also 
an average over p for each of these, for p = 0.1(0.1)0.9 and n = 10, 30, 50. Table 
3 contains 3D approximate theoretical and observed values for Ezs , while Table 4 
contains these values for Ezz , where p = 0.1(0.1)0.9 and n = 10, 30, 50; the second- 
order correction terms for the theoretical values are also stated to 3D. Table 5 con- 
tains 4D values of the observed variance of zs and 3D values of its observed standard 
deviation, likewise for Table 6 with zx , where p = 0.1(0.1)0.9 and n = 10, 30, 50. 
Table 7 contains values of x’ for goodness of fit tests of the normality of zs and zx 
for n = 30, 50. Table 8 contains 2D and 3D values of 


(Ez, — Ez)/+/var (a) + var (2) 


for z, and z, representing the same correlation coefficient but with different p values 
(pz = pi: + 0.1); this is for the product moment correlation coefficient, Spearman’s 
coefficient and Kendall’s coefficient with p, = 0.1(0.1)0.8 and n = 10, 30, 50. 


J. E. Wasa 


Operations Research Group 
System Development Corporation 
Santa Monica, California 


52(K].—G. Horsne.u, “Economical acceptance sampling schemes,” Roy. Stat. 
Soc., Jn., sec. A., v. 120, 1957, p. 148-201. 


This paper is concerned with acceptance sampling plans designed to minimize 
the effective cost of accepted items produced under conditions of normal production. 
Effective cost per accepted item is defined to be the production cost per lot plus the 
average cost of inspection per lot when apportioned equally over the average num- 
ber of items accepted per lot from production of normal quality. Single-sample 
plans are examined in detail. Double-sample plans are considered briefly. 

An appendix contains thirty-one separate tables for single-sampling plans, 











292 REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 


ten of which are applicable when inspection is non-destructive and the remaining 
twenty-one are applicable when inspection is destructive. The table for non-destruc- 
tive inspection displays c,,/c, which is the ratio of manufacturing cost to inspection 
cost; n, the number of items to be sampled; k, the accepted number; A(n, k), the 
probability of acceptance; and c/c,, which is the ratio of effective cost to manu- 
facturing cost. In the tables for destructive inspection, A(n, k) is replaced by 
A'(n, k), which is the expected number of accepted items per lot. 

Plans for non-destructive inspection are given only for a nominal lot size of 
10,000. Plans for destructive inspection are given for lot sizes of 10,000 and 20,000. 
For non-destructive inspection, the process average po = .01(.01).04, the con- 
sumer’s risk point p,; = .03(.01).07, .09; at which the consumer’s risks are .05 and 
.01. For destructive inspection pp = .01 and .02; p; = .03(.01).06, and consumer’s 
risks are .05 and .01. The tables, however, do not include all possible combinations 
of the above listed parameter values. 

A. C. Couen, Jr. 
University of Georgia 
Athens, Georgia 


53(K].—N. L. Jounson, “Optimal sampling for quota fulfillment,’ Biometrika, 
v. 44, 1957, p. 518-523. 


This article contains two tables to assist with the problem of obtaining a preset 
quota m; of individuals from each of k strata by selecting first a sample N of the 
whole population and then completing quotas by sampling from separate strata. 
Individual cost in the first case is c and in the second c; . Table I gives for m; = m 
optimal values of N for k = 2(1)10; mk = 50, 100, 200, 500; d = c;/e = 1.25, 
1.5(.5)3.0; d’ = c//e = .9, .7, .25,0. Here c,’ is the worth of first sample individuals 
in excess of quota. The tabulated values of N are solutions of the equation 
Pr(N; < m) = (ce — e/)/(e: — ¢;). 

Table 2 gives ratio of minimized cost to cost of choosing the whole sample by 
sampling restricted to each stratum. This quantity is 


1 a’ 1\"" (N\ ,, a 
atG-$)0-3) ()e-» 


and is tabulated for k = 2(1)5, 10; km = 50, 100, 500; d = 1.5, 2.5, 3; d’ = .5, 
em 

W. J. Drxon 
University of California 
Los Angeles, California 


54(K|.—P. G. Moors, ‘The two-sample t-test based on range,”’ Biometrika, v. 44, 
1957, p. 482-489. 


This paper provides a sample statistic for unequal sample sizes for a two-sample 
t-test based on observed sample ranges instead of sums of squares. The statistic 
used by the author is simply 
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where Z, and Z, are the sample means, and w, and w, are the sample ranges. Since 
unequal sample sizes are now permitted, the mean range (w, + w.)/2 proposed by 
Lord [1] in his original paper is no longer used. Moore shows that there is very little 
loss in power resulting from the use of the simple sum w, + w:, rather than a 
weighted sum of sample ranges, although for the estimation of o separately, he 


gives a table of f(m , m) and d,, + fd,, to 3D for m , mz. = 2(1)20 to estimate 
o from 


gad wi + fur 
dn, + dn, 


which minimizes the coefficient of variation of the range estimates of population 
standard deviation. 

The main use of this paper is, of course, the tables of percentage points (10%, 
5%, 2%, and 1% points) to 3D for the statistic u, above. The tables of percentage 
points were computed by making use of Patnaik’s chi-approximation for the distri- 
bution of the range, which resulted in sufficient accuracy. The limits for sample 
sizes are m , m2 = 2(1)20, which fulfills the most usual needs in practice. With this 
work of Moore, therefore, the practicing statistician has available a very quick and 
suitably efficient procedure for testing the hypothesis of equal means for two normal 
populations of equal variance. 


F. E. Grupss 
Ballistic Research Laboratory 


Aberdeen Proving Ground, Maryland 


Lorp, ‘‘The use of range in place of standard deviation in the t-test’’, Biometrika, 
v. 34, 1947, p. 41-67. 


55(K}.— NaTIoNAL Bureau or Stranparps, Tables of the Bivariate Normal Distribu- 
tion Function and Related Functions, Applied Mathematics Series, No. 50, 1959, 
xliii + 258 p., 27em. U. S. Government Printing Office, Washington, D. C. 
Price $3.25. 


These tables, compiled and edited by the National Bureau of Standards, pro- 
vide values for the probability content L(h, k, r) of an infinite rectangle with vertex 
at the cut-off point (h, k) under a standardized and centered bivariate distribution 
with correlation coefficient r: 


q pen 1 “ vs oe 2 pe ns 2 A 
L(h, byt) = =a f [ exp | (2 +y aray)/(1 7) |ax ay 


The range of tabulation is h, kK = 0(.1)4, r = 0(.05)0.95(.01)1, the values of 
L(h, k, r) being given to 6 decimal places. For negative correlations, the range of 
tabulation is h, k = 0(.1)hn, kn, r = 0(.05)0.95(.01)1, the values of L(h, k, r) 
being given to 7 decimal places, where L(h,n, kn, —r) < }4.10~ if h, and k, are 
both less than 4. The two tables of L(h, k, r) for positive and negative r, respectively 
(Tables I and II in the text), may therefore be regarded as extensions of Karl 
Pearson’s tables of the bivariate normal distribution in his celebrated Tables for 
Statisticians and Biometricians, Part II, since the range of parameters in the latter 
tables is h, k = 0(.1)2.6, r = —1(.05)1. In this connection, the authors of the 
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present tables present a list of 31 errors in Pearson’s tables, together with the 
corresponding correct values. 

Table III gives the values to 7 decimal places, with the last place uncertain by 2 
units, of the probability content, V(h, \h), of a certain right-angled triangle under a 
centered circular normal distribution with unit variance in any direction. The 
triangle has one vertex at the center of the distribution, with the angle at that 
vertex arc tan A, while the lengths of the two bounding sides at the vertex are h 
and h+/1i + 22. Formally, 


1 Ah Ax r ’ 
V(h, Ah) = xf ax | exp | -5 (x + v) Jay. 


The range of tabulated values for V(h, \h) is h = 0(.01)4(.02)4.6(.1)5.6 and ~, 
and } = 0.1(.1)1. Table IV gives the values of 


2 1 AA afr oy 2 . 
V(Ah, h) = = dx exp 5 (a + y) | dy 


for the parameters h = 0(.01)4(.02)5.6 and ©, and \ = 0.1(.1)1, the degree of 
accuracy being the same as for the values of V(h, \h) in Table III. 

Finally, a short table (Table V) for values of y = are sin r/2x, r = 0(.01)1, 
correct to 8 decimal places is provided. 

The two-parameter function V is related to the three-parameter function L 
by the formula 


L(h,k,r) = V (i, za) +V («, 4%) + F, 





Vi-# Vi-® 
where 

F = 3{l — a(h) — a(k)) + y 
and 


1 4 1 » 
a(x) = —— ex (-5 4) a 
V 29 Jz . : 


This relationship enables L(h, k, r) to be computed to 6-decimal accuracy from the 
7-decimal values of V(h, \h) in regions where interpolation in L(h, k, r) is difficult. 

The function V is also of considerable intrinsic interest, and finds applications 
in such fields as the probability content of polygons when the underlying distribu- 
tion is bivariate normal, the distribution of range for normal samples of size 3, the 
non-central ¢-distribution for odd degrees of freedom, and one-dimensional heat 
flow problems. These applications are illustrated clearly and in detail in the section 
of the book headed Application of the Tables, due to Dr. D. B. Owen. The latter 
section also provides illustrations of the use of L(h, k, r) in problems ranging from 
measurement errors, calibration systems, double-sample test procedures, and per- 
centage changes in sample means to some interesting problems in selection (involv- 
ing the correlation between aptitude test and job performance) and estimation of 
correlation. Finally, an introductory section discusses the mathematical properties 
of the Z- and V-functions as well as methods of interpolation in the tables. 

The authors are to be heartily commended for this most useful book, which will 
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place many statisticians, both practising and those more theoretically inclined, in 
their debt. The visual appearance and general presentation of the material are 
excellent. Perhaps one very minor flaw is that since L(h, k, r) = L(k, h, r), tables of 
L(h, k, r) for h = k would have been sufficient. However, this is a flaw (if at all) 
from the point of view of economics, but hardly so from the point of view of the 
user of the tables! The cost of the book is remarkably low. 


Harotp RvuBEeN 
Columbia University 


New York, New York 


56(K].—E. Nreverce tr, “Die Rangkorrelation U,” Mitteilungsblatt fir Math. Stat., 
v. 9, 1957, p. 196-232. 


In contrast to Spearman’s rank correlation R and Kendall’s coefficient T, the 
author studies the van de Waerden coefficient U. Let p; and q; (i = 1, 2,3 --- n) 
be the ranks of n observations on two variables x and y, and let £;, 9; , ¢; be the 
inverses of the normal probabilities: F(é;) = p:/(n + 1); F(a) = qi/(n + 1); 
F(¢:) = i/(n + 1). Then, U is defined by U = }°2. t/a ¢?. 

If x and y are independent the expectation of U is zero and its standard devi- 
ation is zy = (n — 1)~”, as for Spearman’s coefficient. The author calculates 
also the 4th and 6th moments of U and R, which differ. For n large, U is asymp- 
totically normally distributed about mean zero with standard deviation ov . The 
distribution of U (to 4D) is tabulated completely for n = 4, and over the upper 
5% tail for n = 5, 6, 7. For larger values of n the Gram-Charlier development is 
used. Tables testing independence based on 5%, 2.5%, 1%, and .56% probabilities 
are given to 3D for n = 6(1) 30. For n > 30 the normal probability function 
can be used. 

In the case of dependence the correlation between R and U decreases slowly 
with n increasing. If x and y are normally distributed with zero mean, unit stand- 
ard deviation and correlation p a generalization U* of U to the continuous case 
leads to U* = p. A consistent estimate for p is given. The U test is more powerful 
than the R test. 


E. J. GuMBEL 
Columbia University 
New York, New York 


57(K].—_D. B. Owen & D. T. Monk, Tables of the Normal Probability Integral, 
Sandia Corporation Technical Memorandum 64-57-51, 1957, 58 p., 22 x 28 em. 
Available from the Office of Technical Services, Dept. of Commerce, Washing- 
ton 25, D. C., (Physics (TID-4500, 13th Edn.), Price $.40. 


The following forms of the normal probability integral 
ee a [. edt hzO 
V2 Lx 


G(—h) = 1 — G(h) 


are given for h = 0(.001) 4(.01) 7, to 8D. For those having frequent use of G(h) 
these tables eliminate the simple yet troublesome computation necessary when 
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using the more comprehensive tables prepared at the New York Mathematical 
Tables Project [1] for 





h 
F(h) = OR a 
h 


1 
V 20 
The present tables are more comprehensive than the Pearson and Hartley tables 
[2] for G(h) up to h = 7. These tables have been checked by the authors against 
other tables. (The reviewer could undertake no systematic checking, but such oc- 
casional checks as were made revealed no errors. ) 

Computations (made on a CRC-102A digital computer) were facilitated by 
using the following continued fractions: 

















x h S 25; 
R(x) = To 
34+22 
5— 32? 
. fe 
Rls) 'o 1 ok aa. 
z+1 
r+2 . 
zr+3 
2+: 
Then G(h) was computed from 
” 1 —h2 
G(h) = = 5+ Se PR(h), — ¢"R,(h). 
) Vix 1 - 2 
For large values of h, G(—h) can be obtained from the formula 
G(—h) = eM (h), h=0, 





a 


where 
M(h) =e" [oe *” at. 
h 


M(h), Mill’s ratio [3], provides constants a, an integer, and b (0.1 < b < 1) such 
that 


G(—h) = b.10™ 


For large h = 50(1) 150(5) 500, M(h) and b are tabulated to 8S; a, of course, 
exactly. 

The authors checked —logio G(—h) given in [2] for large h, and found one dis- 
crepancy, namely, —logio G(—500) should read 54289.90830. 


S. B. LirravEer 
Columbia University 
New York, New York 


1. Tables of Normal Probability Functions, National Bureau of he Applied Math. 
Series, No. 23, U. 8S. Government Printing Office, Washington, D. C., 1953. 
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2. E. 8. Pearson, H. O. Harter, aioe Tables for Statisticians, v. 1,.Cambridge 
Universe Press, London, 1954, p. 104-11 


D. Gorpon, “Values of Mill’s atis of area to bounding ordinate of the normal 
probability integral for large values of the argument,”’ Ann. Math. Stat., v. 12, 1941, p. 364-366. 


58(K].—A. E. Sarnan & B. G. Greensere, “Tables for best linear estimates by 
order statistics of parameters of single exponential distributions from singly 
and doubly censored samples,” Amer. Stat. Assn., Jn., v. 52, 1957, p. 58-87. 


Tables are provided for the exact coefficients of the best linear systematic sta- 
tistics for estimating the scale parameter of a one-parameter single exponential 
distribution and the scale and location parameters of a two-parameter single 
exponential distribution. All possible combinations of samples of size n with the 
r, lowest and r, highest values censored are considered for n S 10. Exact coefficients 
for the best linear systematic statistic for estimating the mean (equal to the loca- 
tion parameter plus the scale parameter) are also given for the two parameter case. 
Other tables give the variances, exact or to 7D, of the estimates obtained and the 
efficiency relative to the best linear estimate to 4D based on the complete sample. 
These extensive tables are of immediate practical importance in many fields, such 
as life testing and biological experimentation. 


C. A. BENNETT 
Hanford Laboratories Operation 
General Electric Company 
Richland, Washington 


59(K].—Y.S. Sarne & A. R. Kamar, “Approximations to the distributions of some 
measures of dispersion based on successive differences,’ Biometrika, v. 44, 1957, 
p. 349-359. 


Let 2, --+ , %, be a random sample from a normal population with variance 
o and let 





* 1 n—l 
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The problem is to develop approximations to the distributions of these four types 
of statistics. Let u be any one of these statistics. The method followed is to assume 
that u is approximately distributed as (x,'/c)*, where x, has a chi-square distribu- 
tion with » degrees of freedom; that is, taking \ = 1/a, that cu” is approximately 
distributed as x’ with v degrees of freedom. The constants ¢, a (or \), aud » are then 
determined by equating the first three moments of u to those of (x,’/c)*. The 
results show that a fixed value can be used for @ (or A) if nm 2 5. This allows two 
independent measures of variability u and v2 , based on the same type of statistic, 
to be compared by use of the F test when n 2 5 for both statistics. The basic results 
of the paper are given in Table 1. There, for each of &/o", d/c, 6:'/o°, and d/o, 
fixed values are stated for \, while 3D values for vy and 4D values for logy ¢ are given 
for n = 5(1) 20, 25, 30, 40, 50. Table 2 deals with an example. Table 3 lists the 
results of some approximations to &/o° by (x, /c)* for n = 5, 10, 20, 30, 50. Table 
4 lists for comparison purposes, the upper and lower 1% and 5% points for four 
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approximations to 3'/o° when n = 15, 20. Table 5 is important; it contains 2D 
values of upper and lower 0.5%, 1.0%, 2.5%, and 5% points for the approximate 
distribution developed for ’/o”. Table 6 lists the results of some approximations to 
d/a by (x,’/c)* for n = 5, 10, 20, 30, 50. Finally, Table 7 furnishes 4D values of 
the 6; , 6: differences that result from using a fixed \ for the (x,’/c)* approximation 
to the distribution of 6,’/c’, and from using a fixed for the (x,’/c)* approximation 
to the distribution of d,/c, for n = 5, 7, 10(5)30, 40, 50. 
J. E. Wash 

Operations Research Group 


System Development Corporation 
Santa Monica, California 


60[K].—C. C. Sekar, S. P. Acarwata & P. N. Cuaxrasorty, “On the power 
function of a test of significance for the oo between two proportions,” 
Sankhya, v. 15, 1955, p. 381-390. 


The authors determine the power function of the following statistical test: a 
sample of size n is drawn from each of two binomial distributions with unspecified 
probabilities of success p, and p2, respectively. The null hypothesis is Hy: p, = 
p2 = p. For the two-sided test (alternative hypothesis: p,) < ps2 or pi > pz) at sig- 
nificance level a, the critical region is determined by the following conditions: 

1) For a given total number r of successes in the two samples, the conditional 
probability of rejection under Hy is Sa. 

2) If the partition (a, 7 — a) of r successes is contained in the critical region and 
0 < a <r — a, then the partition (a — 1, r — a + 1) is contained in the critical 
region. 

3) If the partition (a, 7 — a) is contained in the critical region, the partition 
(r — a, a) is contained in the critical region. 

A similar definition is used for the one-sided test of Ho against the alternative 
71 > po. The critical region is determined using the exact conditional probabilities 
for these partitions given by S. Swaroop, [1]. 

The power function for the two-sided test is given to 5D for p; and p, = .1(.1).9; 
nm = 5(5)20(10)50, 100, 200, and for a = .05. For the one-sided test the power 
function to 5D is given for the same levels of p; , po and n, and for a = .025. 

The critical region used by the authors is the one defined for the exact test by 
E. S. Pearson, [2]. However, for small sample sizes the power differs considerably 
from Patnaik’s determinations, which are based on an approximately derived 
critical region and which use a normal distribution approximation of the prob- 
abilities. 

Examples of the use of the tables are included. 


M. L. Epiine 
National Bureau of Standards 


Washington, District of Columbia 


1. Satya Swaroop, ‘‘Tables of the exact values of probabilities for testing the significance 
of pea between proportions based on pairs of small samples,’’ Sankhya, v. 4, 1938, 
p. 7 

2. E. S. Pearson, “The choice of statistical tests illustrated on the interpretation of data 
classed in a 2 x 2 table,’ ’ Biometrika, v. 34, 1947, p. 139-167. 

3. P. B. PATNAIK, “The power function of the test between two proportions in a 2 x 2 
table,’’ Biometrika, v. 35, 1948, p. 157-175. 
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61(K|.—B. Suerman, “Percentiles of the w, statistic,” Ann. Math. Stat., v. 28, 
1957, p. 259-261. 





The statistic 
1 n+1 | 1 n 
.)- ft Bee! Sais EL, 
O36 | *” ati| OSMmS Ti 
is one of several which have been suggested in connection with the null hy pothesis 
that z;, (¢ = 1,---, m) is a random sample from the uniform distribution and 


the L, are the lengths of the n + 1 subintervals of the unit interval defined by the 
ordered sample. In most cases of interest, the z; are the probability transforms of 
observations on a random variable with a continuous distribution function. Based 
on the distribution function derived by the author [1], the 99th, 95th, and 90th 
percentiles of w, to 5D for n = 1(1)20 have been computed, and are given in Table 
I. Values of two standardized forms of this statistic (based on the exact and asymp- 
totic mean and variance, respectively) which are asymptotically normal are given 
in Table II to 5S for the same percentiles as in Table I and for n = 5(5)15(1)20. 
The author points out that the rate of convergence to the limiting values is slow. 


C. A. BENNETT 
Hanford Laboratories Operation 


General Electric Company 
Richland, Washington 


1. B. SnermMan, ‘‘A random variable related to the spacing of sample values,’”’ Ann. Math. 
Stat., v. 21, 1950 p. 339-361. 


62(K, X, Z).—E. D. Casuweitt & C. J. Everert, A Practical Manual on the Monte 
Carlo Method for Random Walk Problems, Pergamon, 1959, 152 p., 23 cm. 
Price $6.00. 


This is volume I of the publisher’s series “International Tracts in Computer 
Science and Technology and Their Application’’. It is devoted to a direct and ele- 
mentary attack on the Monte Carlo principle (that is, the principle of using simula- 
tion for calculation and recording the sample statistics obtained from the simulation) 
in random walk problems, such as mean free path and scatter problems. Many 
examples are given in the text, and an appendix is added listing twenty more or less 
typical problems in which the Monte Carlo method was used at the Los Alamos 
Scientific Laboratory. 

Computational details are given and in many cases flow charts are included. 
No full machine codes are given, but most of the calculations were done on the 
MANIAC I computer at Los Alamos, and coding from the descriptions given and 
the flow charts is probably easier than any attempt to translate a MANIAC I code 
to a code suitable for another machine. A disappointingly short chapter on statisti- 
cal considerations is included; the reader should be warned that this is not a suitable 
exposition of the theory or even the practice of the statistical handling of the sta- 
tistics gathered in his simulation. However, it also is treated from a definitely com- 
putational point of view, including flow charts, and is interesting from this point 
of view. 

An interesting small chapter titled “Remarks on Computation”’ is also included. 
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There is a section on scaling, a section on debugging, a section on special routines, 
a section on a Monte Carlo device for determining the square root of r, and a sec- 
tion on a Monte Carlo device for the cosine of an equi-distributed angle. The ran- 
dom number routine is the familiar routine of selecting the middle digits of the 
square of a quasi-random number; it is frowned on by many random-number special- 
ists. The logarithm routine presented depends on the power series expansion of the 
log, with restriction of the size of the arguments to assure fast convergence. An 
exponential routine is given as a quadratic approximation with scaling of the argu- 
ment. A cosine routine is given through the use of a trigonometric identity and a 
truncated power series for the sine of a related angle. There is no detailed discussion 
of the accuracy of any of these routines. 

The exposition in this book is far from perfect, and the editors have included the 
following statement: ‘It is realized that many workers in this fast moving field 
cannot devote the necessary time to producing a finished monograph. Because of 
their concern for speedy publication, the Editors will not expect the contributions 
to be of a polished literary standard if the originality of the ideas they contain war- 
rant immediate and wide dissemination.” The reviewer feels that the present book 
in its present form is more than justified on the basis of this philosophy, and he 
recommends the book as a most valuable contribution to numerical analysis. 

The chapter headings follow: 

Chapter I. Basic Principles 

Chapter II. The Source Routine 

Chapter III. The Main Free Path and Transmission 

Chapter IV. The Collision or Escape Routine 

Chapter V. The Collision Routine for Neutrons 

Chapter VI. Photon Collisions 

Chapter VII. Direction Parameters After Collision 

Chapter VIII. Terminal Classification 

Chapter IX. Remarks on Computation 

Chapter X. Statistical Considerations 

Appendix. Summary of Certain Problems Run on MANIAC I. 

C. B. T. 


mz 


63(L].—Centre Nationa D’Erupes pes Téiécommunications, Tables numé- 

riques des fonctions associées de Legendre. Fonctions associées de premiére espece, 

P,,” (cos @), deuxiéme fascicule, Editions de la Revue Optique, Paris, 1959, xii + 
640 p., 31 cm. Price 5600 F. 


The first volume of these Tables was reviewed in MTAC, v. 7, p. 178. The 
present second volume was designed to extend the range of tabulation from @ = 90° 
to 6 = 180°. In the process of constructing these tables, however, it was found 
desirable to increase the number of decimals and to add second and fourth central 
differences, thus facilitating interpolation. For this reason, the range up to @ = 90°, 
already covered in the first volume, is included (in an improved form) in the volume 
under review. Perhaps because of the increase in size consequent upon increased 
numbers of decimals and added differences, tabulation has been restricted to m = 
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0, 1, 2 (whereas vol. 1 has m = 0(1)5). Otherwise, the range and interVals of this 
volume match those of the first. 

Since P,,” (cos @) has a singularity at @ = 180° (except when n is an integer), 
an auxiliary function T,,” (cos @) is introduced by the relation 


P,."(cos @) = (ese” @)T,"(cos 0) + (—1)"(Azq logio (cot ‘) P,,"(cos(180° — @)) 
where A, = (2 sin nx)/(x logy ¢). The function 7,” (cos @) is tabulated for 135° < 
6 <= 180°. The use of these auxiliary functions is facilitated by the provision of 
tables of A, , esc 0, esc” 0, and logy cot (6/2). 

The introductory material contains formulas, an account of the tables, hints for 
interpolation, and level curves of P, (cos @), P,' (cos @), and P,. (cos @). 


A. Erpiiyr 
California Institute of Technology 


Pasadena, California 


64(L].—Louis Rosin, Fonctions sphérique de Legendre et fonctions sphéroidales, 
tome 3, Gauthier-Villars, Paris, 1959, viii +- 289 p., 24 cm. Price 5500 F. 


The first two volumes of this work were reviewed in MTAC, v. 13, p. 325f. 
The present, final, volume contains chapters VII to X. 

Chapter VII is devoted to the addition theorems of Legendre functions. Both 
Legendre functions of the first and second kind are included, and two cases are 
distinguished according as the composite argument lies in the complex plane cut 
from — « to +1 or else on the cut between —1 and 1. Addition theorems are also 
developed for the associated Legendre functions of the first kind. 

Chapter VIII is devoted to zeros of Legendre functions. First the zeros of P,.”() 
as functions of yu, for fixed real m and n are discussed, then the zeros of P™,;,(u) 
when m is an integer and p a fixed real number, and then the zeros of Q,”"(). 
This chapter contains also a discussion of zeros of Legendre functions considered as 
functions of n, m and yw being fixed. (These zeros are of importance in certain 
boundary-value problems. ) 

In Chapter [X, applications of Legendre functions are given to partial differen- 
tial equation problems relating to surfaces of revolution other than spheres. Prolate 
and oblate spheroidal harmonics, toroidal harmonics, and conal harmonics are 
discussed. 

Chapter X contains the discussion of some functions related to Bessel functions, 
namely, Gegenbauer polynomials and functions, and spheroidal wave functions. 

Appendix I summarizes relevant information on “spherical Bessel functions’, 
and Appendix II lists numerical tables of Legendre functions and tables connected 
with these functions. 

The third volume maintains the high standards set by the first two volumes, 
and the author must be congratulated upon the completion of this valuable work. 


A. Erp&ty1 
California Institute Technology 
Pasadena, California 
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65(L. M.|.—R. G. Setrrince & J. E. Maxrie.p, A Table of the Incomplete Elliptic 
Integral of the Third Kind, Dover Publications, Inc., 1959. xiv + 805 p., 22 em. 
Price $7.50. 


The advent of fast automatic computers has made a considerable difference to 
the art of making and publishing tables. It has speeded up the computing processes 
fantastically, without, except in relatively minor ways, modifying the labor and 
care needed during publication processes. It has also increased the care needed in 
planning calculations; the plan has to be quite precise and exact in all details for 
the machine to produce proper results, whereas in desk computation, the plan can 
be built up and modified as the work proceeds. 

The glamour of fast computation has led quite a number of people to enter the 
table-making field; people who appear to imagine that the whole problem is simpli- 
fied by automatic computers, who perhaps do not even realize the need to seek 
expert advice. It is with some reluctance, but with the feeling that it is an urgent 
duty that needs to be performed on general grounds, that I suggest that the table 
now reviewed presents one of the most deplorable examples of inadequate planning 
and poor execution that I have met. 

The tables give entries that purport to be 6-decimal values of 





II(¢, a’, k) _ [ 2 2 a ————————————e 
0 (1—a siné)V/1 — Kk sin’? 
for k=sin6, 6=.1(.1)1.5, $ = 0(.01)1.57, 


a = —1(.05) — .1(.02) — .02, .05(.05).5(.02).8(.01).99. 


Also given are two lines for ¢ = 1.5707963, representing a direct and a check calcu- 
lation. 

One example of bad planning is illustrated by the arguments @, and the heading, 
erroneously labelled a. These end arbitrarily in 4 or 5 zeros or 4 or 5 nines. The latter 
is obviously wrong and quite intolerable in print. The present authors are not 
unique in exhibiting this fault, which is quite inexcusable. It is a matter of, at most, 
a few minutes to modify a program, on any automatic machine, to round-off at 
the appropriate figure and to suppress printing thereafter, giving better-looking 
and more convincing argument values. The authors have, in any case given special 
treatment to the argument @ which has one more decimal than function values; 
why not treat it properly? 

Both bad planning and poor execution are exhibited by the check up. The fore- 
word states that “the greatest difficulty was encountered not in constructing the 
table, but in obtaining satisfactory checking’’. This has not, in fact, been achieved! 
The method described could have been—but is not—satisfactory for a” < 0 and 
for 0 < k’ < a’, but the method as described, of integrating through a singularity 
when a’ > 0, is absurdly inadequate. It is not surprising then, that the final two 
lines, both for ¢ = 1.5707963 as mentioned before, should often be in disagreement. 
It is surprising, however, that the authors accept this as a legitimate problem to 
hand on to their readers. The discrepancies indicate errors; the authors’ duty is 
to find and remove these. This states the obvious, but it is equally obviously neces- 
sary to do so. Let it be said again, that to make a program to deal with this properly 
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is a job to be done once and for all. It may take a little effort, time, and money, but 
this is not to be compared with effort, time, and money wasted if it is not done before 
publishing a book of inadequate tables. 

I believe, for my part, the electronic computation is so fast and easy that a dis- 
crepancy of more than a single unit (known or unknown) is not tolerable in a 
published table (of which the computing effort and cost are now often only a small 
part of the whole effort and cost); however, I am willing to concede that discrep- 
ancies up to perhaps two units might, on rare occasions, be justified, provided this 
is clearly stated. It is quite intolerable to have the following discrepancies; to pick 
out some of the worst: 














Page a2 k2 l(a’, &) 
| 

661 0.87 0.86869... 8.654098 and 8.654259 
733 0.93 0.92844... 15.309882 and 15.310251 
745 0.94 0.92844... 16.857725 and 16.857859 
793 0.98 0.97111... 45.498015 and 45.498457 
805 0.99 0.92844... 49.243943 and 49.244046 
805 0.99 0.97111... 70.018520 and 70.018897 





The discrepancies are, in fact, highly systematic throughout; they are all of the 
same sign, except for the violent cases mostly listed above near the singularity 
k = a mentioned in the introduction. They indicate clearly that at least one set 
of the check values is erroneous because of an inadequate method, and not merely 
because of rounding; severe doubt is cast on both discrepant values. Dr. J. W. 
Wrench, Jr. has computed anew the value on p. 733 for a’ = .93, k’ = sin’ 1.3, and 
finds 15.3098662, which is not even between the two values quoted from the book; 
neither published value is correct and one errs by more than the discrepancy be- 
tween them. I repeat again, it is the duty of the table compiler to remove all these 
doubts. 

The poor execution is also exhibited by the fact that in the heading, a and k 
appear in place of the correct a’ and k’, while the 10-decimal values of k’ given 
(which are simply sin’ \ for \ = .1(.1)1.5) have end figure errors running up to 11 
units. Again, the argument is given as @ in the tables; this corresponds to ¢ in the 
introduction. It is only fair to add that the heading errors in a’, k’, @ have been 
announced as errata. Another awkward point for the user is that absence of values 
for a = 0, and for k = 0, makes the tables harder to interpolate. 

From all this, it is evident that the authors are lacking in experience of table- 
making so that their remark ‘With the argument as outlined, no attempt has been 
made to proof or check the printed sheets in any way other than by comparison of 
the resultant complete integrals” causes less surprise than might otherwise be the 
case. If there were no other faults occurring other than those mentioned above, 
they would have been exceptionally and unduly lucky. Photographic processes 
seem as far from infallibility as printing from letter-press; the possible faults are 
different, but nevertheless exist just the same. In fact, a rather superficial exami- 
nation of the table reveals unsightly irregularities in spacing of lines on pages 51, 
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623, 700 and some lesser ones elsewhere. It would have been easy to reprint the 
pages before reproduction. More serious are several digits that are not fully legible: 


p. 446 6 = .25 k = .03946 --- 3rd digit 3 
p. 507 6 = .75 k = .15164 --- 1st digit 8 
p. 579 Bottom right corner. Very faint. 


All these imperfections occur in at least two copies of the tables; such imperfections 
are common in tables printed from typescript and should be expected and sought 
out. The real surprise is, however, as mentioned above, that, after the numerical 
comparison of check values mentioned had been made, its lack of success seems 
simply to have been ignored. 

It is hoped that possible users may, with the exercise of necessary—but undue— 
caution, obtain adequate results, maybe 5} correct figures, if they need them. The 
publication of this book will undoubtedly make it much more difficult to publish a 
good and proper version; this is a major criticism of such a book. The only consola- 
tion I can offer the authors is that I have seen several tables that are even worse. 

As I have said, I have expressed myself so freely with some reluctance, from a 
sense of duty; it is no part of my desire to discourage the enthusiasm of table-makers, 
but they must realize the magnitude and duties of the task so taken on, and seek 
competent advice before proceeding with the work, and potential users must be 
adequately warned. 

J. C. P. Minuer 
The University Mathematical Laboratory 
Cambridge, England; and 
The Mathematics Research Center, U. 8. Army 
University of Wisconsin, Madison, Wisconsin. 


66(S].—D. R. Hartree, The Calculation of Atomic Structures, John Wiley and Sons, 
Inc., New York, 1957, xiii + 181 p., 23 cm. Price $5.00. 


This book by the late D. R. Hartree is the fruit of a lifetime of experience in the 
calculation of the outer, electronic structure of atoms. It is concerned with methods 
for the calculation of atomic structures rather than with the results of such calcula- 
tions for particular atoms. Emphasis is deliberately placed on means of obtaining 
“best” approximations which can be both represented and applied simply. The 
student who wants an introduction to the essential methods of approximation and 
computation of shell structures may read the first hundred pages. The mathemati- 
cian will find in this book the physical background for the author’s well-known text 
on numerical analysis. 

In the Introduction are outlined the seven main steps in the development of 
atomic theory up to the point at which quantitative calculations are possible. The 
atomic units are introduced and the point change approximation of the electron 
justified. Then, properties of the Schroedinger equation are summarized to prepare 
the reader for the main problem of the book, the numerical solution of the self- 
consistent field equations with and without exchange. The variation principle is 
carefully introduced, and the total energy of closed shell configurations discussed. 
Also, configurations comprising incomplete groups are treated. In the later part, 
the main ideas and methods are extended to more complicated or more complete 
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cases, and, in general, are described only very briefly. The text concludes with a 
chapter on “Better Approximations”. 

Many tables, relating to Slater coefficients, mean radii, screening numbers and 
reduced radial wave functions, are found in the text and in Appendix 2. For the 
bibliography up to 1947, the author refers to Reports on Progress in Physics, v. 11, 
1946-47, p. 141-143, and completes the list to October 1956 in Appendices 1 and 3. 

Numerical procedures, most of them recommended by the author’s experience 
in hand and machine calculations, are described in detail, giving step-by-step 
instructions and numerical examples. The reviewer would have liked to have seen 
some comments on the numerical stability of methods using finite differences 
approximations to differential equations. Numerical stability is obvious in the 
Numerov and Fox-Goodwin process, but this is not so in Hartree’s method of para- 
graph 4.6 (page 72), although the application is correct. Familiarity with numerical 
stability prevents the physicist from blindly refining the methods given in the text 
and will save costly “numerical experimentation”’. 

In general, the book would gain by stating briefly the mathematical reasons why 
certain procedures are recommended (it would be mostly in the light of numerical 
stability!), as, for example, for the separation of integration of the radial wave 
equation into an outward and inward integration (§5.2). The physical reasons are 
stated adequately. 

It seems that many equations were renumbered before the manuscript went to 
the printer. Cross references to equations are quite often unreliable. Otherwise, the 
number of misprints for a book of this kind is rather low. 


Erwin BareEIss 
Argonne National Laboratory 
Lemont, Illinois 


67(W, X].—Rosert O. Fereuson & Lauren F. Sarcent, Linear Programming, 
McGraw-Hill Book Co., New York 1958, xiv + 342 p., 24 em. Price $10.00. 


With increasing frequency the professional mathematician, especially if he is 
working in applied mathematics, finds himself approached by friends or colleagues 
who lack advanced mathematical training but want to know more about the latest 
techniques, such as linear programming. In the course of the past year, several 
books on linear programming have appeared to which such inquirers might be 
referred. 

This book should probably be regarded as the best of the group. It is addressed 
primarily to “people engaged in management activities at all levels in the firm and 
students of management... .” Its major virtues include a simple expository style 
without condescension, a wealth of illustrative examples, and a somewhat broader 
coverage of the subject than other works currently available. As a result of these 
qualities, it should prove suitable for individual study by management personnel 
with substantial practical experience in industry. It should, however, have its 
greatest value as a textbook for classroom instruction (on the job or off) of groups 
in which some or all participants lack the mathematical prerequisites which would 
permit use of a more advanced text, such as the well-known volume by Gass (from 
the same publishing house, interestingly enough). 

The three sections into which the book is divided are entitled Introduction, 
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Methods, and Application, while two technical appendices treat the mathematical 
foundations of the Simplex Algorithm, and its relationships to the modi method, 
respectively. The second section presents the transportation method, the modi 
method, the simplex method, all of which are exact, and two approximation 
methods: the index method and the authors’ own ratio-analysis method. Applica- 
tions illustrated in detail in the third section include a product-mix problem, a pro- 
duction smoothing problem, and a problem in optimal assignment of orders to 
plants, where production costs as well as distribution costs affect the decision. 
Emphasis in this section is on obtaining and utilizing a maximal amount of informa- 
tion of value to management from the linear programming solution. Other applica- 
tions, it may be added, are illustrated in presenting the computational methods in 
the second section. Further, two technical appendices treat the mathematical 
foundations of the Simplex Algorithm and its relationships to the modi method, 
respectively. 

The professional mathematician should probably be warned that, while he may 
safely recommend this book to non-mathematicians, he should not attempt to read 
it himself, unless he is willing to take the risk of apoplexy. As evidence for this con- 
clusion, one may cite a number of quite apoplectic reviews in which this book and 
others like it have been severely castigated (by mathematicians) for a low level of 
scholarship, lack of rigor, and similar mortal sins. To some, this attitude may seem 
unfair. Texts on business statistics, to take an analogy at random, are not usually 
criticized for omitting discussion of Borel sets, Stieltjes integrals, and the Cramer- 
Rao inequality. One might expect that there would be a place for a comparable 
treatment of linear programming without reference to theorems on matrix inversion 
or convex polyhedral cones. 

For good or for ill, linear programming is being dished up for the common folk, 
and this book represents probably the most workmanlike presentation currently 
available. As might be expected, the book is weakest where it is most technical. The 
attempt, on page 50, to explain degeneracy explains nothing. Experience indicates 
that a better treatment of this concept is possible without resorting to advanced 
mathematics. Similarly, on page 5 and again on page 77, the difference between 
linear programming and the solution of simultaneous equations is explained in 
terms of the non-optimizing character of the latter, but the authors do not go on to 
explain, as they easily might have done, the reasons for this difference. Also, one of 
the examples (p. 119 ff.) used to illustrate the simplex method contains one 
redundant equation (because it is inherently a transportation problem) but the text 
makes no mention of this fact. As is all too often the case in technical works, the 
index is far from complete. For example, under ‘“‘degeneracy”’ there is no reference 
to page 50, which has the only non-technical discussion of the topic (such as it is). 

All these omissions could easily be remedied, and the many good qualities of the 
book warrant the hope that a future edition will see such improvements made. 


Tuomas A. GOLDMAN 
Corporation for Economic and Industrial Research 
Arlington, Virginia 
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68[X].—Kurt Arbenz, Integralgleichungen fir einige Randwertprobleme fitr Gebiete 
mit Ecken. Promotionsarbeit. Eidgendssische Technische Hochschule, Zirich, 
1958, 43 p. 


This paper is devoted to the problem of finding a procedure, suitable in numeri- 
cal application, for the conformal representation of a simply connected plane domain 
over the unity circle, in the more difficult case of a boundary with corners. A useful 
tool which has been employed by Todd [1] is the integral equation of Lichtenstein 
and Gerschgorin, but it cannot be directly applied in this case. Modifications have 
been proposed by Stiefel and by Birkhoff, Young and Zarantonello [2]. 

The author generalizes theorems given by Radon, on the potential theory for 
domains bounded by smooth ares. He makes extensive use of methods of functional 
analysis, with particular reference to the book of Riesz and Nagy [3]. 

The last seven pages contain numerical examples: the conformal representation 
of a square on the unit circle, obtained in a rather simple way with good accuracy; 
and the displacements in a square plate with built-in boundary. It appears that no 
use has been made of electronic computers, and it would be of interest to start 
numerical experiments on computers with the procedure here suggested. 


Enzo L. Aparo 
University of Rome 
Rome, Italy 


1. NBS Apritiep MatuemarTics Series, No. 42, Experiments in the Computation of Con- 
formal Maps, U. S. Government Printing Office, Washington, D. C., 1955. 

2. Proceedings of Symposia in Applied Mathematics, v. 4, 1953, p. 117-140. 

3. F. Rresz & B. V. Sz. Naay, Legons d’analyse fonctionelle, Third Edition, Gauthier- 
Villars, Paris, 1955. 











TABLE ERRATA 


283.—G. W. Spencetey, R. M. Spencetzey & E. R. Epperson, Smithsonian 
Logarithmic Tables to Base e and Base 10, Smithsonian Institution, Washington, 
D. C., 1952. 


page entry for read 
2 In 86 4.45435.--- 4.45434-.-- 
39 In 1931 7 .56475--- 7 .56579- - - 
221 log 915 95142-.-- 96142--- 
288 log 4271 63052 94714--- 63052 95714--- 


C. R. Sexton 
2947 Elmwood Court 


Berkeley 5, California 


Eprtori1aL Nore: For additional errata, see MT'AC, v. 10, 1956, p. 261; v. 11, 1957, p. 226. 


J. W. W. 


284.—E. Campi, Eleven and Fifteen-Place Tables of Bessel Functions of the First 
Kind, to All Significant Orders, Dover Publications, New York, 1948. 
On p. 27, the value given for Jo(9.76) should read —.22836--- instead of 
— .23836---. 


C. R. Sexton 





285.—W. MEYER zuR CaPELLEN, [ntegraliafeln. Sammlung unbestimmter Integrale 
elementarer Funktionen, Springer, Berlin, 1950. [See RMT 1090, MT AC, v. 7, 
1953, p. 100.] 

In addition to errata listed on a loose sheet supplied with the book and those 
appearing in MTE 223 (MTAC, v. 7, 1953, p. 106), the following corrections 
should be made: 

p. 245, column 3, formula 1.5.1.1, line 1. 





for 008 §2 [Ci{k(a + 2)} + Ci{k(a — 2))}* 
read or [Ci{k(a + x)} — Ci{k(a — 2)}]*. 





p. 245, column 3, formula 1.5.1.2, line 1. 


for — 998 Ea ici{k(e + a)} + Ci{ (2 — a)}}* 
reed == [Ci{k(a + a)} — Cifk(a — a)}]*. 


Mur.an 8S. Corrineton 
Radio Corporation of America 
Camden, New Jersey 
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NOTES 
New Journals 

Two new journals of interest to applied mathematicians have recently been 
announced. 

Journal of Mathematical Analysis and Applications, published by Academic 
Press, New York and London. This new journal will provide a medium for the 
rapid publication of carefully selected mathematical papers treating classical 
analysis and its manifold applications. 

The refereeing system used in most journals is replaced by a board of Associate 
Editors, each of whom may accept manuscripts. In this manner the delay between 
receipt and publication of manuscripts will be minimized. Only papers written in a 
lucid, expository style will be eligible for publication. 

In recognition of the fact that other disciplines contribute new concepts and 
problems to the continuing growth of mathematics, papers devoted to the mathe- 
matical treatment of questions arising in physics, chemistry, biology, and enginéer- 
ing will be encouraged. In these papers the emphasis will be upon the analytical 
aspects and the novelty of problem and solution. 

Journal of Mathematical Physics, published by the American Institute of Physics, 
New York, is a bimonthly devoted to new mathematical methods for the solution 
of physical problems as well as to original research in physics furthered by such 
methods. Its scope includes: mathematical aspects of quantum field theory, statis- 
tical mechanics of interacting particles, new approaches to eigenvalue and scatter- 
ing problems, theory of stochastic processes, novel variational methods, theory of 
graphs, and review papers on mathematical topics for physicists. 


Automatic Programming of Digital Computers—National Information Centre 


The National Centre of Information on Automatic Programming of Digital 
Computers has been established by the Department of Mathematics of the Brighton 
Technical College, Brighton, England, in response to a recommendation of the 
first National Conference on Automatic Programming, held in Brighton in April, 
1959. This conference was attended by 111 delegates from computer manufac- 
turers, industrial and commercial computer users, government research establish- 
ments, universities and technical colleges. 

The purpose of the Centre is: 

(i) to establish and maintain a comprehensive library of publications, papers, 
research reports and other material, especially those not readily accessible, in 
any way relevant to the problems of automatic programming, and to make these 
available, in English, on demand; 

(ii) to publish, in conjunction with Pergamon Press Limited, an annual review 
in Automatic Programming: 

(iii) to provide a ‘clearing house’ for information and inquiries on automatic 
programming and related topics, and to help co-ordinate the work of other bodies 
active in this field; 

(iv) to organize from time to time, small working conferences on particular 
aspects of the subject; 

(v) to maintain permanent contact with organizations in all other countries 
concerned with such matters. 
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CORRIGENDUM 


2 2 P —¢2 
R. Hensman & D. P. Junxins, Tables of =e" i e~" dt for Complex z, Math. 
Comp., Review 12, v. 14, 1960, p. 83. 


for read 


2 22 [ —t2 2) 2% [ —t2 
of. 5% 7 A o dt 
The research was done at the Royal Radar Establishment, Malvern, Worcester- 
shire, England. 








San, 


